,^J^^^' 


.V\o?,' 


Courant  Institute 

of  Mathematical  Sciences 

MAGNETO-FLUID  DYNAMICS  DIVISION 


MP-36  NYO-10^3^ 

Finite  Amplitude  Compression  Waves 
In  a  Collision-Free  Plasma 

K.W.  MORTON 
15  January  1964 
AEC  Research  and  Development  Report 


\ 


NEW    YORK    UNIVERSITY 


-»:,»,»^ 


4»    #^  l^- 


This  report  was  prepared  as  an  account  of  Government  sponsored  work.  Neither 
the  United  States,  nor  the  Commission,  nor  any  person  acting  on  behalf  of  the 
Commission: 

A.  Makes  any  warranty  or  representation,  express  or  implied,  with  respect  to 
the  accuracy,  completeness,  or  usefulness  of  the  information  contained  in 
this  report,  or  that  the  use  of  any  information,  apparatus,  method,  or 
process  disclosed  in  this  report  may  not  infringe  privately  owned  rights;  or 

8.  Assumes  any  liabilities  with  respect  to  the  use  of,  or  for  damages  resulting 
from  the  use  of  any  information,  apparatus,  method,  or  process  disclosed 
in  this  report. 


As  used  in  the  above,  "person  acting  on  behalf  of  the  Commission"  includes 
any  employee  or  contractor  of  the  Commission,  or  employee  of  such  contractor, 
to  the  extent  that  such  employee  or  contractor  of  the  Commission,  or  employee 
of  such  contractor  prepares,  disseminates,  or  provides  access  to,  any  information 
pursuant  to  his  employment  or  contract  with  the  Commission,  or  his  employment 
with  such  contractor. 


UNCLASSIFIED 


Magneto-Fluid  Dynamics  Division 
Courant  Institute  of  Mathematical  Sciences 
New  York  University 


TID-4500 

26th  Edition  Physics 


MF-56  NYO-10^3'4- 

Finite  Amplitude  Compression  Waves 
in  a  Collision-Free  Plasma 

K.W.  MORTON 
15  January  196'!- 
AEC  Research  and  Development  Report 


Contract  AT (30-1 ) -1^80 


-  1  - 
UNCLASSIFIED 


ABSTRACT 

Adiabatlc  two-fluid  equations  are  set  up  to  describe 
the  one -dimensional  motion  of  a  collision-free  plasma.   These 
are  then  used  in  a  study  of  finite  amplitude  compression 
waves  which  is  built  around  the  consideration  of  a  "magnetic 
piston  problem" . 

The  asymptotic  behaviour  of  small  amplitude  waves  is 
considered  first.   It  is  found  to  give  a  simple  but  qualita- 
tively correct  picture  of  how  the  flows  depend  upon  the  two 

parameters  9,    the  angle  the  initial  magnetic  field  makes  with 

2 
the  direction  of  wave  propagation,  and  e  ,  the  electron-ion 

mass  ratio.   When  tan  0  >  £  -  e,    the  main  wave  front  is 

followed  by  an  oscillatory  wave  train  on  the  scale  of  the 

geometric  mean  of  the  ion  and  electron  gyromagnetic  radii; 

otherwise,  the  front  is  preceded  by  an  oscillatory  precursor 

on  the  scale  of  the  ion  gyromagnetic  radius. 

An  exhaustive  study  is  carried  out  of  the  steady  flows 

which  are  useful  in  the  solution  of  the  piston  problem.   Their 

2 
character  depends  markedly  on  the  parameters  9,    e   and  the 

Alfven  Mach  number  M. .   The  relationship  to  the  corresponding 

one-fluid  flows  emerges  and  the  effect  of  introducing  friction 

forces  is  also  considered. 

Finally,  the  piston  problem  is  solved  numerically  using 

finite  differences.   The  results  clearly  show  that  the  behaviour 

deduced  from  the  asymptotic  theory  and  the  steady  flow  is 

substantially  correct,  even  quantitatively  so  in  many  cases. 
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1 .   Introd-uctlon  and  Summary 

The  study  of  shocks  and  finite  amplitude  waves  in  a 
collision-free  plasma  has  received  recent  stimulus  from  several 
quarters.   For  example,  one  proposed  method  of  heating  a  plasma 
is  to  pass  shocks  through  it.   But  in  a  hot  rarefied  plasma  the 
mean  free  path  may  be  so  large  that  the  ordinary  collisional 
shock  becomes  a  rather  ineffective  means  of  heating.   In  order 
to  pursue  this  idea  then,  one  has  to  seek  a  shock  mechanism 
Independent  of  collisions.   In  another  field,  the  possibility 
that  waves  of  this  type  might  exist  in  the  solar  wind  has  been 
suggested.   However,  once  these  problems  are  raised,  they  have 
sufficient  inherent  interest  to  maintain  their  own  momentum, 
and  we  shall  not  discuss  the  motivation  for  their  study  any 
further. 

One  of  the  fundamental  problems  concerns  the  existence  of 
a  shock  in  a  two-component,  collision-free  plasma  containing  a 
magnetic  field.   Even  if  one  confines  oneself  to  essentially 
one-dimensional  phenomena,  this  involves  the  self -consistent  solu- 
tion of  the  Vlasov  equations  for  the  two  types  of  particle  together 
with  Maxwell's  equations.   Some  progress  in  this  direction  has  been 
made  by  Morawetz  [25].    Specifically,  she  shows  that  the  shock 
would  have  an  oscillatory  fine  structure  and  proposes  as  the 
governing  mechanism  the  "bouncing"  of  the  ions  against  the  result- 
ing potential  barriers.   But  the  parameter  ranges  in  which  her 


See  also  Molseev  and  Sagdeev  [23]. 
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treatment  is  valid  are  narrow  and  a  simpler  model  Is  needed 
for  a  broad  survey  as  well  as  for  the  solution  of  practical 
problems. 

Thus  one  turns  to  macroscopic,  i.e.   fluid,  descriptions 
of  the  plasma.   The  simplest  such  model  is  a  one-fluid 
model  in  which  the  plasma  is  regarded  as  a  perfectly  con- 
ducting fluid  with  an  adiabatic  pressure.   The  resulting 
equations  were  derived  by  Lundquist  [21]  and  have  yielded 
a  very  rich  theory  which  has  been  extensively  studied  (see 
Grad  [15]  for  a  survey).   However,  the  model  gives  no  shock 
structure  at  all  and  has  little  obvious  relevance  for  a 
collision-free  plasma. 

Between  these  two  extremes  lie  many  models,  one  of  which 
is  studied  in  the  present  work.   This  is  the  adiabatic  two- 
fluid  model  proposed  and  briefly  discussed  by  Gardner  et  alia 
[11].   The  total  mass,  momentum  and  energy  conservation  equa- 
tions are  as  in  the  one-fluid  model  but  the  individual 
momentum  equations  for  the  two  types  of  particle  are  combined 
to  give  an  equation  for  the  current;  this  eliminates  the  need 
for  assuming  any  kind  of  Ohm's  law.   The  model  has  the 
advantage  of  still  being  quite  a  simple  fluid  model  but  also, 
in  the  important  case  of  a  cold  plasma,  of  being  completely 
equivalent  to  a  self -consistent  particle  formulation.   Indeed 
it  serves  to  relate  several  disparate  elements  of  collision- 
free  shock  theory  (see  Blank  and  Grad  [?]  on  this  subject). 

One  cannot,  of  course,  expect  a  proper  shock  structure 


from  such  a  model  for  It  Includes  no  specific  dissipative 
mechanism.   But  one  can  obtain  non-trivial  compression  waves 
of  finite  amplitude,  and  for  a  considerable  range  of  amplitudes 
these  do  not  break  to  form  shocks.   Even  when  shocks  do  form 
one  can  hope  to  localise  them  and  hence  those  regions  in 
which  a  more  exact  theory  is  necessary.   As  in  classical  gas 
dynamics  the  canonical  problem  is  that  of  the  compression 
wave  generated  by  the  movement  of  some  sort  of  piston.   It  is 
around  such  a  problem  that  our  study  will  be  built. 

As  a  cold  plasma  theory,  the  two-fluid  model  has  been 
studied  by  several  authors,  especially  the  transverse  case, 
in  which  the  magnetic  field  is  everywhere  normal  to  the  direc- 
tion of  wave  propagation.   Attention  was  focussed  on  steady 
one-dimensional  flows,  and  for  this  special  case  a  solitary 
wave  and  periodic  solutions  were  bound  by  Adlam  and  Allen  [l], 
Davis,  Liist  and  Schlttter  [9]^  and  others.   More  recently, 
Saffman  considered  first  the  longitudinal  case  [32]  and  then 
the  general  case  of  arbitrary  initial  orientation  of  the  mag- 
netic field  [33]-   The  initial  state  constitutes  a  singular 
point  of  the  equations  and  he  concentrated  his  attention  on 
the  nature  of  this  point,  finding  the  conditions  under  which  a 
solution  may  grow  from  it.   At  about  the  same  time,  first 
Morawetz  [24]  and  then  Ban'os  and  Vernon  [4]  and  Vernon  [35] 
considered  the  full  equations  with  the  pressure  included  and 
studied  steady  solutions  in  the  transverse  case,  but  they  still 
restricted  themselves  to  cases  in  which  no  breaking  occurred. 
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The  work  of  these  authors  Is  limited  by  Its  considera- 
tion of  only  steady  flows.   The  fundamental  problem  Is  the 
time -dependent  piston  problem  and  the  steady  flows  should 
arise  naturally  from  this;  the  same  basic  principle  is 
embodied  in  the  evolutionarity  condition  of  Polovin  [29]  for 
the  existence  of  shocks  in  conservative  systems.   By  failing 
to  adopt  such  a  point  of  view,  one  is  likely  either  to  obtain 
unstable  steady  flows  or (as  in  the  present  case)  to  miss  a 
significant  proportion  of  the  meaningful  solutions.   But  to 
solve  the  present  time -dependent  problem  analytically  is  out 
of  the  question;  the  most  that  has  been  done  in  this  direction 
is  the  asymptotic  solution  for  weak  non-linear  waves  by 
Gardner  and  Morikawa  [13].   So  the  answer  lies  in  using 
numerical  solutions  of  the  time -dependent  problem  to  guide 
one's  search  for  relevant  steady  flows.   This  is  the  method  of 
attack  used  here;  the  results  for  the  transverse  case  were 
published  In  an  earlier  paper  (Morton  [26])  and  the  present 
work  covers  the  case  of  arbitrary  initial  field  orientation. 
Its  success  is  such  that  sufficient  steady  flows  are  found  to 
describe  nearly  all  parameter  choices  for  the  piston  problem 
and  the  numerical  solutions  can  almost  be  dispensed  with.   Auer, 
Hurwitz  and  Kilb  [2],  [3]  have  also  studied  the  time -dependent 
problem  for  the  transverse  case  numerically  but  with  a  semi- 
physical  model  which  aimed  at  obtaining  the  shock  structure. 
They  did  not  attempt  to  link  their  results  with  the  steady  flows 
of  the  two-fluid  model. 


The  plasma  Is  assumed  to  be  initially  uniform,  at  rest,  ' 
and  containing  a  uniform  non-zero  magnetic  field.   The  flows 
to  be  studied  are  one  dimensional,  directed  along  the  x-axis. 
At  the  beginning  of  section  2  the  equations  of  motion  are  set 
up  and  made  dimenslonless .   As  a  result  the  problem  is  found 
to  depend  on  the  following  five  parameters: 

6,  the  angle  between  the  initial  magnetic  field  and  the 

X-axis , 

2 
8  ,  the  ratio  of  the  electron  to  the  ion  mass, 

M.,  the  Alfven  Mach  number, 

P,  the  ratio  of  the  fluid  to  the  magnetic  pressure 

initially,  and 

7,  the  adiabatic  exponent  of  the  total  fluid  pressure. 
In  most  of  the  work  (all  the  numerical  work)  we  have  taken 

Y  =  2,  the  appropriate  value  for  the  transverse  case.   Also, 
since  we  know  from  the  work  of  Gardner  and  Lynn  [12]  and  others 
that  the  adiabatic  pressure  is  not  a  very  good  approximation 
when  the  plasma  is  warm,  we  have  usually  assumed  that  it  is 
initially  cold,  and  so  taken  p  =  0.   It  may  be  remarked  here 
that  our  object  in  including  any  pressure  terms  at  all  is  not 
to  give  an  accurate  description  of  a  warm  plasma,  but  rather 
to  allow  a  correct  qualitative  picture  of  flows  after  breaking 
has  occurred.   Hence  we  take  the  simplest  scalar  pressure; 
various  anisotropic  pressure  tensors  have  been  used  by  Rose 
[31]  J  Hain,  LUist  and  Schltiter  [l6],  and  others,  but  the  results 
have  shown  the  same  major  features  in  the  cases  studied. 
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After  a  discussion  of  the  general  properties  of  the 
equations,  small  amplitude  waves  are  studied  In  ^2.^.      The 
dispersion  relation  for  the  linearised  equations  Is  derived 
and  then  an  asymptotic  analysis  like  that  of  Gardner  and 
Morlkawa  [13]  Is  carried  out.   The  result  Is  quite  striking. 
It  shows  that,  when  tan  0  >  e~   -  e,  weak  compression  waves 
have  an  exponentially  rising  front  followed  by  an  oscillatory 
wave  train;  while  for  tan  0  <  e    -  e   an  oscillatory  precur- 
sor runs  ahead  of  the  main  front  which  now  rises  smoothly  to 
a  final  constant  state.   Further,  the  wave  lengths  of  the  two 
sets  of  oscillations  are  very  different;  that  of  the  precur- 
sor scales  like  the  Ion  gyro-radius,  while  the  wave  train  Is 
on  a  scale  sometimes  called  the  "magnetic  penetration  depth" 
but  which  Is  actually  equal  to  the  geometric  mean  of  the  Ion 
and  electron  gyro-radll  (the  corresponding  frequency  Is  also 
sometimes  called  the  hybrid  frequency).   These  features  prove 
to  anticipate  the  whole  body  of  results. 

In  the  next  and  longest  section,  section  5?  steady  flows 
are  studied.   The  steady  state  equations  In  general  reduce  to 
two  second  order  differential  equations  In  the  two  transverse 
components  of  the  magnetic  field.   As  they  are  written  In 
equations  (5-9)  to  (3.12),  they  Involve  the  specific  volume  V 
which  In  turn  Is  given  by  a  quadratic  (3-13)  whose  coefficients 
are  functions  of  the  fields  and  their  first  derivatives.   The 
roots  of  this  quadratic  determine  two  sheets  on  which  solutions 
may  lie.   And  the  regions  In  phase  space  In  which  these  roots 
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are  real  and  positive  play  a  vital  role  In  limiting  the 
existence  of  physically  meaningful  solutions.   Shocks  occur 
as  solutions  jump  from  one  sheet  to  the  other. 

The  singular  points  of  the  steady  state  equations  are 
the  initial  state  and  all  the  states  behind  the  expansion 
waves  and  the  fast,  slow,  and  intermediate  shocks  obtained 
from  the  de  Hoffman -Teller  shock  relations.   These  and  the 
solutions  in  their  neighbourhood  form  the  subject  of  *5  3-2. 
It  is  followed  by  a  study  of  the  steady  flows  in  three  special 
cases  in  which  a  particularly  detailed  treatment  is  possible 
because  of  reduction   of  the  equations.   In  two  of  them,  the 
transverse  case  and  the  equal  mass  case,  the  equations  for 
the  y-  and  z-components  of  the  fields  dissociate.   The  solu- 
tions are  similar  in  the  two  cases,  being  of  the  wave  train 
type,  and  those  for  the  former  case  are  described  in  the  earlier 
paper,  Morton  [26].   The  third  case,  in  which  the  electron 
mass  is  neglected,  does  not  seem  to  have  been  studied  before. 
Here  the  two  steady  state  equations  reduce  to  first  order. 
The  existence  of  a  solitary  wave  beginning  and  ending  at  the 
fast  shock  singularity  and  circling  the  initial  state  is 
proved.   Inside  are  the  periodic  waves  which  give  rise  to  the 
precursor  type  solutions  usual  in  this  case.   In  9J>A   the 
general  case  is  shown  to  be  intermediate  in  behaviour  between 
these  extremes;  the  transition  is  described  and  is  seen  to  be 
centered  roughly  on  the  line  tan  9   =   e~      -  e.   Finally,  the 
effect  of  introducing  a  friction  term  into  the  equations  is 
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considered.   It  Is  shown  that  In  some  cases  the  initial  and 
shock  states  can  then  be  joined  by  a  continuous  solution. 

The  piston  problem  is  considered  in  the  fourth  and  final 
section.   The  motion  of  the  piston  is  here  simulated  by  the 
turning  on  of  either  an  electric  or  magnetic  field  at  a 
fixed  plane  wall  marking  the  boundary  of  the  plasma.   The 
field  is  raised  gradually  to  a  final  value  in  order  to 
initiate  a  steady  flow  in  the  plasma.   Finite  difference 
equations,  using  an  artificial  viscosity  to  allow  accurate 
approximation  to  shocks,  are  used  to  solve  this  time -dependent 
problem.   The  anticipated  steepening  of  the  compression  waves 
takes  place  until  an  almost  steady  wave  form  is  attained.   In 
most  cases  this  can  be  matched  quite  accurately  with  steady 
flows  determined  by  the  Initial  state  and  fast  shock  singu- 
larities; as  might  be  expected  the  accuracy  is  worst  for  weak 
waves  and  for  parameter  values  in  the  transition  zone  near 
tan  9   =   E        -   e. 

The  main  results  are  for  p  =  0.   They  show  that  for  Mach 

numbers  up  to  a  value  which  ranges  from  )/2~  to  2  according 

2 
to  the  values  of  e   and  0,    there  is  no  breaking  of  the  wave 

2 
and  no  completely  steady  flow.   Starting  from  e  =  1  or 

9   =   7r/2  the  main  feature  of  solutions  is  a  wave  train  on  the 

scale  length  of  the  mean  gyro-radius  already  described.   As 

the  transition  region  is  reached  the  precursor  grows  while 

the  wave  train  decreases  in  amplitude  as  in  Figure  12.   Then 

2 
at  the  other  extremes,  0  =  0  or  e  =0,  the  precursor,  scaling 
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like  the  Ion  gyro-radius,  Is  the  dominant  feature.   Higher 
Mach  numbers  bring  a  change  In  the  form  of  the  solution  behind 
the  main  front.   The  chief  effect  occurs  when  M.  =  Mp(0),  a 
function  which  varies  from  Y^  at   9   =   0   to  2  Aj   at  0  =  7r/2 ; 
above  this  value  the  fast  shock  and  the  Initial  state  singu- 
larities are  on  different  V-sheets  and  the  solutions  must  con- 
tain a  discontinuity  In  V.   Further,  at  these  Mach  numbers 
the  nature  of  the  fast  shock  singularity  has  changed  so  that 
at  one  extreme  the  wave  train  disappears  while  at  the  other 
a  wave  train  Is  added  to  the  precursor.   In  some  of  these 
cases  genuinely  steady  flows  are  approached  as  t  ^^oo. 

These  results  show  that  even  for  a  cold  plasma,  physically 
realistic  values  of  the  mass  ratio  result  In  the  Ion  gyro- 
radius  being  the  dominant  length  as  soon  as  one  departs  a 
little  from  the  transverse  case.   Similarly  Gardner  and 
Morlkawa  [1^],  working  with  the  Vlasov  equation,  found  that 

the  weak  solitary  wave  In  the  transverse  case  was  broadened 

1   ~2   1/2 
by  the  factor  (l  +  y^  e  p )  '^   when  the  temperature  was 

Increased.   Together  the  two  results  suggest  that  the  cold 
transverse  case,  which  for  a  long  time  received  most  of  the 
attention.  Is  very  special.   As  soon  as  Its  strong  restric- 
tions on  the  motion  of  the  electrons  are  relaxed,  by 
Increasing  the  temperature  or  adding  a  longitudinal  component 
to  the  magnetic  field,  the  oscillations  at  the  hybrid  frequency 
are  damped  out  and  the  Ion  gyro-radius  becomes  the  length  scale. 
The  relationship  between  the  one-  and  two-fluid  theories 
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is  also  established.   The  states  ahead  and  behind  each  shock 
of  the  former  theory  correspond  to  singular  points  of  the 
equations  In  the  latter.   In  some  cases  these  can  be  joined 
by  a  steady  flow;  then  the  relationship  Is  clear.   But  alsoj 
the  parameter  ranges  where  this  Is  not  possible.  I.e.,  when 
the  oscillatory  wave  train  or  precursor  Is  present,  shrink  as 
P  Increases.   Indeed  as  f3  — ^  os  all  solutions  become  of  one 
type  In  which  the  fluid  variables  behave  as  In  the  one-fluid 
theory  while  the  magnetic  field  has  a  continuous  monotonia 
profile  as  In  Figure  7- 

Thus  the  jump  In  the  magnetic  field  which  occurs  In  all 
the  one-fluid  shocks  is  resolved  by  the  two-fluid  theory. 
The  most  Important  question  which  remains  unanswered  Is  how  to 
resolve  the  discontinuities  In  the  fluid  variables  which  still 
sometimes  occur.   For  this  theory  to  have  any  value  this  must 
be  accomplished  by  a  mechanism  operating  on  a  length  scale  at 
most  no  greater  than  those  of  the  two-fluid  solutions.   We  do 
not  attempt  to  answer  this  question,  but  two  points  are  worth 
noting.   First,  Zwlck  and  Gonzales  [56]  have  computed  self- 
consistent  particle  orbits  In  the  case  e   =1,0  =  7r/2  up  to 
very  high  Mach  numbers.   The  orbits  are  looped  and  If  a 
mechanism,  such  as  the  two  stream  instability,  were  able  to 
smooth  out  these  loops  it  appears  that  there  would  be  a  marked 
similarity  to  the  solutions  of  Figure  6.   Secondly,  work  done 
by  Liubimov  [20]  and  Hu  [private  communication]  indicates  that 
it  is  incorrect  to  assume  that  in  the  presence  of  a  magnetic 
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field  a  small  number  of  collisions  would  lead  to  a  shock  with 
a  width  of  the  order  of  the  mean  free  path.   In  particular, 
the  former  author  obtains  a  width  of  the  order  of  the  ion  gyro- 
radius  in  the  case  0=0. 
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2 .   The  Equations  of  Motion 
2.1.   Derivation 

The  fundamental  equations  governing  the  dynamics  of  a 
collision-free,  two-component  plasma  consist  of  the  Vlasov 
equations  for  the  distribution  functions  f_|.(r,w,t)  of  the 
two  types  of  particle  and  Maxwell's  equations  for  the  electro- 
magnetic fields  E(r,t)  and  B(r,t).   However,  In  studying  only 
low  frequency  fluid  phenomena  we  may  omit  the  displacement 
current  In  the  latter  giving  Instead  the  pre-Maxwell  equations; 
thus  we  have,  writing  S,  for  S/Bt, 

(2.1)      a^f^  +  w.  V/+  +  (e+/m+)[(E  +  w  x  B)  •  V^U]    -  0 

V  X  B  =  M- J 
'2.2,    ^ 

V  ^  E  +  a^B  -  0 

Here  the  suffices  +  and  -  distinguish  between  the  positively 
charged  Ions  and  the  negative  electrons  with  charges  e+  and 
masses  m+;  the  pre-Maxwell  equations  are  In  rationalised 
m.k.s.  units  with  [x   the  permeability  and  k  the  dielectric 
constant  of  free  space.   The  charge  density  q  and  current  J 
are  given  by  the  following  two  moments  of  the  distribution 
functions : 


;2.3)      q  =  <e>  J  =  <ew>; 
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here  the  notation  <(Q,^,  defined  for  any  scalar  or  vector  functlc 
of  e^,  m^.,  and  w.  Is  used  to  denote 


<Q>  =    (Q+f+  +  Q_f_)d^v 


The  two-fluid  equations  on  which  the  present  work  Is 
based  may  be  derived  directly  from  the  above  equations  by 
taking  moments.   The  crucial  approximation,  which  together 
with  the  assumption  of  one -dimensional  geometry  enables  this 
to  be  carried  through,  is  the  neglect  of  the  displacement 
current.   For  then  (2.2)  Implies  that  the  x-component  J-,  Is 
zero,  which  In  turn,  through  the  charge  conservation  equation 
h.q   +  S  J  =0,  Implies  that  If  q  =  0  Initially  It  remains  so. 
Polsson's  equation  Is  then  not  used;  taking  q  =  0  Instead  Is 
an  assumption  usually  called  "quasl-charge  neutrality". 

However,  we  wish  to  use  only  the  simplest  fluid  pressure 
and  equations  of  state,  namely  a  scalar  pressure  p  with  the 
equation  of  state  appropriate  for  a  polytroplc  gas  with 
adlabatlc  exponent  y.   In  this  case  It  Is  a  simple  matter  to 
write  down  the  mass,  momentum,  and  energy  conservation  equa- 
tions directly;  these  lead  to  the  equations  (2.8)-(2.12)  below 
for  the  density  p  (or  specific  volume  V)  and  mean  velocity  Uj 
where 

(2.^)  p  =  1/V  =  <m>   ,      pu  =  <^mw^  . 

The  equations  relating  J  to  the  fields  and  taking  the  place 
of  an  Ohm's  law  are  more  difficult  to  obtain;  since  they  con- 
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stltute  the  essential  difference  of  this  model  from  the  one- 
fluid  model  of,  say,  the  Lundqulst  equations  we  shall  give 
their  derivation.   We  multiply  the  equations  (2.1)  by  ew. 
Integrate  over  all  w-space,  and  add  the  equations  together 
to  get  after  a  little  rearrangement 

(2.5)      S^J  +  a^  <ew^w>  =  (E  +  u  X  B)  (^^   _  b  K  (^  (w-u)^ 


We  need  only  y-  and  z-components  of  this  equation  (to  be 
denoted  by  suffices  2  and  5);  and  to  evaluate  the  averages  we 
have  to  use  the  assumptions  that  q  =  J-,=  0  and  that  the  stress 
tensors  are  diagonal,  which  was  implicit  in  the  assumption  of 
a  scalar  pressure.   We  get 


<(ew^w^>  =  u^J^   ,     1  =  2,3 


(2.6)   {  <e^/m>  =  a^e^ 


P 
<(e^/m)(w-u)>  =  -a(l-£^)j 


where  a  =  -  e_/m  is  the  charge  to  mass  ratio  of  the  electron 
and  e   =  -  m_e,/ni  e_  is  the  electron  to  ion  mass  ratio  in  the 
case  of  singly  ionised  ions.   After  differentiating  (2.5)  with 
respect  to  x  we  can  eliminate  E  and  J  by  using  (2.2).   Then  to 
obtain  the  full  set  of  equations  in  the  form  given  below  we 
change  to  dimenslonless  variables  and  Lagrangian  co-ordinates. 
Velocities  are  scaled  by  the  initial  Alfven  speed  a  =  B  /(m-P 
and  lengths  by  one  of  the  two  alternatives  given  by 


o^ 


;2.7)       { 

2  - 

1/x^  =  ae  (m-Pq)^. 


Here  p   and  B  are  the  initial  density  and  magnetic  field; 

"^  o        O  o  o  ' 

X  may  be  regarded  as  either  the  ratio  of  the  speed  of  light 

the  gyromagnetic  radii  of  the  ions  and  electrons  moving  with 
speed  a  across  the  initial  field;  while  x.  is  the  gyromagnetic 
radius  of  the  ions.   It  will  be  convenient  sometimes  to  use 
one,  sometimes  the  other.   Finally,  normalising  to  a  unit' 
initial  density  and  magnetic  field  and  changing  to  the 
Lagrangian  co-ordinate   J  p dx  which  we  continue  to  denote 
by  X  we  have , 

(2.8)  S^V  -  a^u^  =  0 

(2.9)  ^^{p-v"^)    =   0 

(2.10)  S^u^  +  a^(p  +  I  B^  +  I  b|)  =  0 

(2.11)  a^Ug  -  B^  b^B^    =   0 

(2.12)  b.u^    -  B,  a  B^  =  0 
^    ^         t  3    1  X  3 

(2.13)  B^  -   constant,  cos  9 

(2.14)  h^{B^Y    -   A  a^B^)  -  B-^(a^U2  +  f^^B^)  =    0 


Lagrangian  co-ordinates  are  more  convenient  both  for  applying 
boundary  conditions  and  for  numerical  integration  of  these 
equations.  _  ..  „ 


(2.15)  \(^3^  -  ^  ^x^^)  -  ^l^^x''^  "  r'^A^  =  °- 

The  constants  A  and  \~'   depend  on  the  scale  length  chosen: 

mean  scale  length  x:-     A==l,P'=e    -e 
ion  scale  length  x.:~  A=e,P=l-e 

The  electric  field  E,  normalised  by  a  B  ,  and  the  current 
7,  normalised  by  aea  p  ,  are  given  by  the  subsidiary  equations 

(2.16)  Eg  =  u^B^  -  B^u^  -  Aa^S^B^  +  TB^a^Bg 

(2.17)  E^  =  B^Ug  -  u^Bg  +  AS^a^Bg  +  Pb^^^B^ 


;2.18)  Jg  =  -A^p  a^B^ 


;2.19)  J3  =  A^p  S^Bg 


2.2.   General  Properties  of  the  Equations 

That  the  variables  in  these  equations  (2.8)-(2.12)  are 
far  more  loosely  coupled  than  those  in  the  Lundquist  equa- 
tions soon  shows  itself  in  a  consideration  of  their 
characteristics.   The  fatter  equations  yield  three  signifi- 
cant characteristic  speeds  and  studies  of  the  characteristics 
have  produced  many  important  and  interesting  results.   But 
the  only  non-degenerate  characteristic  speed  of  the  two-fluid 
equations  is  the  adiabatic  sound  speed.   The  remaining 
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characteristics  consist  of  the  particle  paths,  t  =  constant, 
five  times  and  x  =  constant  four  times.   The  last  represent 
the  11ml o  of  two  pairs  of  characteristics  travelling  with  the 
speed  of  light,  as  this  tends  to  Infinity;  this  may  be  seen 
by  using  the  full  Maxwell's  equations  but  retaining  the  assump- 
tion of  charge  neutrality  in  the  derivation  of  the  equations. 
Since  also  we  are  particularly  interested  in  the  case  of  an 
initially  cold  plasma  the  charac;teristics  provide  very  little 
Information;  the  most  that  can  be  said  is  that  they  indicate 
the  correct  boundary  conditions  to  be  applied. 

In  the  Lundquist  equations,  (2.5)  is  replaced  by 
E  +  u  X  B  =  0,  which  results  in  the  reduction  of  (2.l4)  and  (2.I5) 
to 
(2.20)         ^^{BV)    =   B3_B^u  , 

for  the  y,z  components.   This  is  the  only  difference  between 
the  two  sets  of  equations.   Now  it  is  well  known  (see 
Llubarskii  and  Polovln  [I9]  and  Shercliff  [34])  that  compres- 
sive simple  waves  in  the  Lundquist  model  steepen  to  form 
shocks,  so  when  this  is  a  good  approximation  to  the  two-fluid 
model  the  same  phenomenon  should  occur  there.   This  will  be 
the  case  when  the  terms  in  S  B  occurring  in  (2.l4)  and  (2.I5) 
are  small  compared  with  the  others;  how  small  these  higher 

derivatives  have  to  be  will  in  general  depend  on  the  parameters 

2 
e  ,  6      and  particularly  on  the  speed  or  amplitude  of  the  wave. 

Thus  such  a  "gentle"  perturbation  of  an  initial  state  should 

steepen  where  there  is  compression  until  the  diffusive  effects 
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of  the  second  derivatives  take  over.   We  shall  see  below 

that  this  does  indeed  happen  in  the  piston  problem  if  the 

piston  is  not  pushed  too  hard. 

If  the  piston  is  pushed  harder,  however,  shocks  do 

form  though  they  will  differ  from  those  in  the  one -fluid 

models.   For  in  weak  solutions  of  the  equations  we  can  see 

that  for  A  /^  0,  B  and  5  B  are  always  continuous  while  if 

A  =  0,  B  is  continuous.   Thus,  if  for  each  fixed  value  of  x 

we  denote  by  L  the  operator  /  dt ,  then  from  (2.11)  and 
^  o 

(2.12),  we  obtain  for  the  y-  and  z-components 


'2.21)     u  =  B.L  a  B 
'  1  X  X 


Operating  with  L  on  (2.l4)  and  (2.15)  we  get 

(2.22)  (A  +  B^L^)a^B2  +   B-^  Tl^S^B-^     =B2V 

(2.23)  -BirVx^2    +  (^  +  ^l^x^^x^3  =  ^^    -    ^^"  ^ 

valid  even  when  the  right  sides  are  discontinuous.   When  A  /^  0 
this  is  a  coupled  pair  of  Volterra  integral  equations  of  the 

second  kind  for  ^^B^^.   Their  solution  is  simple  and  we  obtain 

2 
the  boundedness  of  S^B  from  that  of  V  which  will  be  assured 

when  there  is  no  cavitation.   In  the  valuable  limiting  case 

A  =  0,  obtained  by  using  the  ion  scale  length  and  letting  the 
2 


ma 


ss  ratio  e   tend  to  zero,  we  need  B-,  /^   0.   Then  the  bounded- 
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2 

ness  of  V  implies  that  of  L  5  B,  so  for  discontinuities 

travelling  at  a  finite  speed  SB  is  bounded. 


2 .3  Small  Amplitude  Waves 

In  this  section  we  consider  only  a  cold  plasma,  and, 
mainly  for  later  reference,  we  first  record  the  dispersion 
relation  for  waves  linearised  about  the  initial  state  V  =  1, 
u  =  0,  B  =  (cos  9,  0,  sin  9).      After  the  elimination  of  all 
other  variables,  the  linearised  equations  for  the  magnetic 
field  become 

(S^  -  cos^e  a^  -  AS^S^)B^  =  pcos  6   S.a^B 


t  x^3 


;2.2^) 


(^'t  -  ^x  -  ^^t^x)^^  =  -  1"^^°^  '    \^l^2 


Substituting  the  plane  wave  B(x,t)  =  B  exp  <  ik(x  +  Ut)( 
yields  the  dispersion  relation 


2^.2  i-i  2    2, 


(2.25)     [U  (1  +  Ak  )-  1][U  (1  +  Ak  )-cos^0]  -  k  U  P  cos 


Note  that  when  9  /=   7r/2  and  the  mean  scale  length  is  used,  the 

2 
behaviour  of  the  wave  velocity  as  the  mass  ratio  e  —^   0  is 

given  by 


'2.26)    U  --  (1  +  k^)  ^k  cos  0  •  e  -^ 
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This  Is  In  agreement  with  the  behaviour  of  the  solitary- 
wave  found  by  Saffman  [32]  In  the  longitudinal  case,  0=0. 

Now  we  turn  to  the  consideration  of  the  asymptotic 
form  of  a  wave  Initiated  by  the  slow  compressive  motion  of 
a  piston.   This  has  been  studied  In  detail  by  Gardner  and 
Morlkawa  [13]  In  the  transverse  case;  we  first  repeat  their 
treatment  for  the  general  case  and  then  return  to  the  matter 
of  Its  justification.   Thus  we  make  perturbation  expansions 
of  the  variables  V,  u-,  ,  u^^  B^  In  the  form 


(0)  ,  R  (1)    r2  (2 
=w^  ^+5w^  ^+5w^ 


+  .  . 


and  of  the  variables  u^,    B^,   which  are  Initially  zero.  In  the 
form 

w  =  52(0  +  Sw^^''  +  S^w^^-*  +  ...). 

Now  It  Is  known  (see  below)  from  the  study  of  the  steady  state 
equations  that  the  speed  of  propagation  of  a  small  amplitude 
compression  wave  Is  approximately  the  magneto-sonic  speed, 
which  Is  unity  here.   So  we  change  to  the  following  "stretched' 
variables, 

(2.27)         I  =  5*(x-t)  ,     T  =  5^/2  t  . 
Equations  (2.8)-(2.13)  then  become 

(5^^-  a^)V  =  S^u^ 

(Sa^-  a^)u^  +  B^'6^B^   +   B^a^B-^  =:  0 
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(65^-   5^)u2   =   cos    e    S^B^ 


5S    -    hi-  )u^   -   cos   0    SftB^ 


2^    N  .       ^/^  .     r-o^^S, 


(6B^-    a^)(B2V    -    SAB^Bg)    =    cos    e{h^u^   +    f  b^hf^B^ 


2^      .  ^/-^  ^o?n2. 


6S^-    S^)(B^V    -    6Aa^B^)    =   cos    0(a^u^    -    P^^^^p, 


'2 


To   lowest    order     u\°^    =  u^°^    =  u^^-*    =  bI^-*    =   0,   V^^-*    =   1, 
I  d  ^  d 

B4    '    =   sin   0;    and   to   first    order  we    obtain 

a^(V^^)   +  u['^h    =   0  ^^(4'^''  -  sin  0B^^))  =  O 

^^^^2       +  cos  0  u^^^  +  rcos0  ^pB^^h  =  0,    Sfi  (u^-^^+  cos  0B^^^)  =  0 
B|(sln0  V^-^^+ B^-^''+  cos  0u^^^)=  0  B^  (u^"^  ^  +  cos  0  B^^  ^ )    =   0 

For   0  <   0  <   ir/2 ,    these   equations  have   the    solutions 


uj^-*    =    -V^"^''    =   sln0B^"^'*    =    -tanOu^^-*    =  y,      say. 


and 


■^^    =    -cos  0B^^''    =    pcos^0    sin   ^0    S^y 


U^     '     =     -cos    -xjp 


Then  a  non-linear  equation  for  y  results  from  requiring  the 
consistency  of  the  2nd  order  equations: 

(2)  (2)  2 

SfiCuj  ^-   sin  0  B^  ^)    =   hY+   cosec  Ojh^y 
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sin  9    5s (sin  9  Y^^'   +   B^^^  +  cos  9   u^^ 


3    '     ---    -    -3 


cos  e  a^y  +  2jb^Y   +  (A  -  r  cot  9)b^Y 


L(u^^)  +  cos  0)  B^^^)  =  -  cot  9    a^y 


Multiplying  these  equations  by  sin  9,    -sin  9,    -cosec  9,    and 
cos  0  respectively  and  adding  we  obtain 

(2.28)         25^y  +  3yS|y  +  (A  -  r^cot^0)S^y  =  0. 


This  is  the  same  equation  as  that  obtained  by  Gardner 
and  Morikawa  [13]  apart  from  the  coefficient  (A  -  ["*  cot  0). 
Let  us  denote  a  typical  solution  of  their  equation  by  16(^,1). 
Then  when 

tan  0  >  £   -  e 

the  corresponding  solution  of  (2.28)  is  given  by  y(l,T)  = 
j6{^    jt')  under  the  simple  scale  change  ^  =  (A  -  |   cot  0)^^  , 
T  =  (a  -  I   cot  0)^t'.   Clearly  the  appropriate  scale  length 
in  this  case  is  the  mean  ion-electron  gyro-radius  in  which 
A  =  1.   Then  y  has  the  same  general  form  as  (6   which  typically 
consists  of  a  progressing  wave  front, rising  exponentially 

from  the  initial  state,  followed  by  a  decaying  wave  train 

1/3 
which  spreads  with  time  like  t  ^     .      On  the  other  hand,  when 

£   -  £  >  tan  0  the  solution  of  (2.28)  is  given  by 

y(^,T)  =  -j?5(-|  '  ,t' ),  where  I  =  (  r^cot^0  -  A)  ^4  '  and 
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'1     'J  -^ 

T  =  ( P  cot  9    -   A)^T  .   Thus  the  wave  front  Is  now  preceded 
by  an  oscillatory  precursor  and  converges  smoothly  behind 
to  the  state  at  the  piston.   Further,  the  change  of  scale 
Is  approximately  proportional  to  "  so  in  this  case,  which 
includes  the  usual  situations  where  e  is  small,  the  appro- 
priate scale  length  is  the  ion  gyro-radius  in  which  |   =  1-e  , 

It  remains  to  give  some  motivation  for  the  initial 
expansions  and  the  change  of  variables  (2.27).   We  shall  see 
in  the  next  section  on  steady  flows  that  a  progressing  wave 
in  which  B^  =  sin  6   +  0(5)  travels  at  a  speed  M»  =  1  +  0(6). 
Moreover,  in  the  neighbourhood  of  the  initial  state 
B^  ^  (B^-  sin  0)(M^  -  l)^  =  0(5-^/^).   Hence  the  different 
expansion  for  Bp  and  Up  is  suggested.   Also,  from  momentum 
and  energy  conservation  for  a  cold  plasma  (see  equations  (3.2 
to  (3-8))  one  obtains 


[h^Bf   =   M^^[(M^  -cos^9)(B2  +  (B^-sin  ef) 


+  B^  -  sin^e)^]  =  0(5^; 


The  appropriate  stretching  for  x  follows,  and  for  the  t-stretchin^ 
we  rely  on  the  results  of  Gardner  and  Morikawa  [13]-   By  taking 
Laplace  transforms  they  obtained  an  asymptotic  solution  of  the 
linearized  equations  (2.24)  for  9   -   it/2   which  depended  solely 
on  the  variable  (x-t)t  '     . 

We  have  given  this  treatment  in  some  detail  because  the 
results  of  later  sections  show  that  it  gives  an  essentially 
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correct  picture  of  how  the  solutions  to  the  piston  problem 
depend  on  the  two  parameters  9   and  e  when  the  disturbance 
Is  weak. 
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3-   steady  Flows 

3.1.   The  Steady  State  Equations 

A  great  deal  of  information,  both  of  general  Interest  and 
of  particular  value  for  the  piston  problem,  can  be  obtained 
from  the  steady  flows  which  satisfy  equations  (2.8)  to  (2.I5), 
These  are  solutions  which  appear  time -Independent  to  an 
observer  moving  with  a  fixed  velocity  M„,  the  Alfven  Mach 
number.   The  ordinary  differential  equations  which  they  satisfy 
are  obtained  by  replacing  time  derivatives  by  -M.5  and  are 
usually  called  the  steady  state  equations  or  equations  of 
steady  motion.   They  can  immediately  be  Integrated  once,  but 
in  order  to  include  shocks  in  this  formulation  we  have  to  be 
careful  to  use  equations  in  their  conservation  form.   Thus 
equation  (2.9)  should  be  replaced  by  the  energy  conservation 
equation 

(3.1)     \e^^^   +  S^(u^p)  =  A-*V  E.J 

where  e,  , ,  the  total  energy  per  unit  mass,  is  given  by 


:3.2)     e^^^  =    (^_1)-V  +  | 


and  V  =  A~^VJ. 

The  constants  of  integration  are  given  by  the  initial 
state,  which  is  also  the  state  at  x  =  +00,  namely. 
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V=l,  ^  =  0,  p  =  Jp  sin  9,  ..     ,.  ..  r:.,i.   .,■" 

(3.5) 

B  =  (cos  6,    0,    sin  0),   E  =  J  =  0     ' .',   -■  '-J    ■■  ■■   -j 

where  p  Is  the  ratio  of  the  Initial  fluid  and  magnetic  pres- 
sures.  Then  one  obtains  from  the  conservation  equations 

(3.4)  u^  =  M^(l-V) 

(3.5)  p  +  I  (B^  +b|)  =  M^u^  +  I  (l+p)sln^9 

(3.6)  M^u^  +  B^B2  =  0 

(3.7)  M^u^  +  B^B^  =  B^sln  0 

(3.8)  "^A^tot  ^^l^^i"  \^^2 +  ^3^  "  V3  ^^"  ^  "  constant. 

To  obtain  the  last  equation  we  use  the  facts  that  Ep  =  M»  (B-^-sln  9) 

and  E-,  =  -M.Bp  which  are  deduclble  from  equations  (2.l6)  and 

(2.17). 

The  first  four  of  these  equations  together  with  the 

definition  of  e,  ,  can  be  used  to  eliminate  p  and  u  from  the 
tot 

last,  thus  yielding  a  quadratic  for  V  with  coefficients  depending 
on  B  and  v.   This  and  the  Integrated  forms  of  (2.l4)  and  (2.I5) 
constitute  the  complete  set  of  steady  state  equations  for  B, 
V ,    and  V : 

dB2 


3.9) 

dx 

^3 

3.10) 

dB^ 

dx     ~ 

-V2 

30  - 


dv       Tb. 


(3-13)     F(B,v)V^  +  G(B,v)V  +  H(B,v)  =  0 


where 


P  =  (7  +  1)M^ 


(3.1^)     G  -  -7[2M^  +  (l+p)sln^e  -  (B^  +  b|)] 

H  =  (y-1)[M^-  2  sin  0(B  -sin  0)  -  B^M~^  (B^  +(B^-sln0)^; 
-  a{v^   +  v^)]  +  yp  sin  9    . 


If  one  considers  only  solutions  leading  from  the  Initial 
state  and  If,  further,  one  assumes  this  state  Is  cold.  I.e., 
^   -   0,    and  that  the  Mach  number  Is  sufficiently  small  for  the 
solution  to  remain  continuous,  then  one  can  work  with  a  simpler 
equation  than  (3-13).   For  then  p  =  0  so  {^.^)    and  (3-5)  give 


(3.15)    V*  =  1  -  (Bg  +  b|  -  sln^0)/2M^ 


This  Is  the  equation  used  by  many  earlier  authors,  e.g., 
Adlam  and  Allen  [l]  and  Saffman  [33]-   To  see  the  relation 
between  the  two  equations  we  substitute  (3.15)  Into  (3'11) 
and  (3-12)  which  will  allow  us  to  carry  out  a  further  Integra- 
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2     2 
tlon.   This  yields  an  expression  for  Vp  +  v^,  which  in  the 

special  case  of  equal  ion  and  electron  masses  (see  below, 

3. 3(a))  virtually  completes  the  determination  of  the  solution. 

In  general,  by  substituting  this  expression  into  (3-1^)  it 

enables  us  to  check  that  V   is  indeed  a  root  of  the  quadratic 

(^•13)-   It  is  also  easy  to  see  that  the  other  root  is 

V  (7-l)/(7+l).   As  we  shall  see  below,  a  switch  of  roots  for 

V  corresponds  to  a  shock;  so  we  have  shown  that  a  shock  from 
a  cold  initial  state  leads  to  the  maximum  compression 
(7+1)/(y-1)  as  in  ideal  gas  dynamics  (see  e.g.,  Courant  and 
Friedrichs  [ 8] ) . 

However,  to  derive  the  maximum  amount  of  useful  informa- 
tion from  a  study  of  the  steady  state  equations  it  is  necessary 
to  adopt  a  more  general  point  of  view.   We  need  to  consider 
solutions  other  than  those  which  arise  directly  from  the 
initial  state,  and  we  wish  to  study  the  occurrence  of  shocks. 
Thus  we  shall  use  the  system  (5-9)-(3-12)  together  with  the 
subsidiary  definitions  (3. 13)  and  (/^.l^).      These  equations  form 
an  autonomous  system  of  ordinary  differential  equations  in  the 
four  dimensional  phase  space  of  the  variables  (Bp,  B^,  Vp ,  v^). 
Rather,  they  form  two  such  systems,  one  for  each  root  of  (3-13) 
for  V;  correspondingly,  we  may  talk  of  solutions  lying  on  one  of 
two  sheets.   The  right  hand  sides  of  the  equations  are  functions 
which  are  everywhere  continuous,  and  indeed  are  analytic  except 
where  the  two  roots  of  (3-13)  are  equal,  i.e.,  on  the  surface 
where  the  two  sheets  intersect.   There  are  mathematical 
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difficulties  In  continuing  solutions  across  this  surface. 
But  since  It  also  forms  the  boundary  of  a  region  In  which  the 
roots  for  V  are  complex,  physical  reasons  alone  are  sufficient 
for  one  to  consider  only  those  solutions  which  remain  In  Its 
complement . 

The  excluded  region  Is  connected  and  Its  form  can  be 

easily  deduced  from  the  expressions  for  F,  G,  and  H.   These 

2       2     2 

depend  on  v  only  through  ^      =   A(vp  +  v-^)  so  that  the  geometry 

Is  essentially  three  dimensional.   The  equation  G  =  0  defines 

a  cylinder  U.  with  axis  parallel  to  the  ^-axls  while  H  =  0 

defines  a  spheroid.   When  projected  onto  a  plane  i^  =  const. 

these  give  two  circles   and  the  V-roots  are  complex  In  an 

annular  or  crescent -shaped  region  centred   roughly  on  or 

Inside  the  former  and  lying  Inside  the  latter;  typical  con- 

2 
figurations  are  shown  In  Figure  1.   As  (^   Is  Increased,  H 

Is  reduced  and  the  projection  of  the  excluded  region  shrinks 

to  zero.   It  Is  worth  remarking  that  for  some  parameter  values 

there  Is  no  excluded  region.   For  example,  this  Is  so  for 

0   =   0,    1  <  M^  <   3/2   and   for  M^  =   1 ,    0  <   0  ^  tt/8. 

Physical  conditions  impose  two  further  restrictions  on  the 

solutions.   For  the  specific  volume  to  be  positive,  solutions 

on  the  sheet  corresponding  to  the  smaller  value  of  V  must  be 

confined  to  the  interior  of  the  surface,  H  ==  0.   And  on  both 

sheets,  solutions  must  be  inside  the  cylinder  G  =  0  if  the 

pressure  is  to  be  positive,  since  it  is  easy  to  see  that 

G  =  -2y(p  +  m|v). 


In  the  transverse  case,  H  =  0  defines  a  straight  line  parallel 
to  the  Bp-axis  for  fixed  (^ . 
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3-2 .   singular  Points  of  the  Steady  State  Equations 
(a)   The  location   of  the  singular  points 

The  first  step  In  studying  the  qualitative  character  of 
the  solutions  of  the  equations  (5- 9) -(5- 12)  Is  to  locate  the 
singular  points  and  determine  their  nature.   They  are  obtained 
by  equating  all  derivatives  to  zero  so  thay  correspond  to 
constant  states  of  the  physical  system  where  variables  will 
be  distinguished  by  a  tllda;  we  have 


(5.16)         v^  =  v-j 


(5.17)  (M^V  -  B^)B^  =  (M^  -  B^)sln 

(3.18)  (M^V  -  B^)B2  =  0 


and  these  equations  together  with  {'^  A)-{J>  .8)    completely  des- 
cribe such  states.   This  set  of  equations  Is  now  Independent 

2 
of  the  mass  ratio  e  and  Indeed  of  this  particular  model,  for 

they  are  equivalent  to  the  de  Hoffman-Teller  hydromagnetlc 

shock  relations  which  have  been  studied  exhaustively  by  many 

authors  Including  Heifer  [17].  Frledrlchs  [10],  Liist  [22], 

Bazer  and  Erlcson  [6],  Shercllff  [34],  and  Polovln  [29]. 

The  constants  have  been  chosen  so  that  one  of  the  singular 

points  always  corresponds  to  the  Initial  state  (3.3).   The 

remainder  correspond  to  final  states  of  various  types  of  waves 

usually  categorised  in  the  above  references  as  expansion  waves. 
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Figure  I.      Typical  Behavior  of  the  Roots  for  the   Specific   Volume 
V  in  the   Plane  7=0- 
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slow  shocks :,  and  fast  shocks  with  Alfven  (or  Intermediate) 
shocks  and  contact  discontinuities  occurring  as  limiting  or 
transitional  cases.   The  characterisation  of  these  waves  is 
quite  complicated  (see  e.g.  Bazer  and  Ericson  [6])  but  can 
be  presented  fairly  simply  for  the  present  purposes.   We 
begin  by  specialising  to  p  =  0  and  7=2.   From  (3-l8),  we 
see  that  in  all  but  certain  special  cases  Bp  =  0;  and  (5-17) 
gives  an  expression  for  V  which  we  substitute  into  the 
quadratic  (3-13) •   The  result  is  a  relationship  between  M. 
and  B^  with  0   as  the  sole  parameter.   Explicitly, 


;3.i9)   M^  =  (3-n)~"^[Ti-3(T]-i)  cos^e  +  n(i-^  {n-ifs±n^2e)^] 


where 


T]  =  B^/sln  9 


3 

and  the  initial  state  has  been  factored  out.   These  roots  are 
shown  in  Figure  2;  for  each  value  of  B   there  is  a  continuous 
looped  curve  which  begins  at  M.  =  0,  B^  =  sin  0,  V  =  1  as  a 
slow  shock,  changes  at  M.  =  cos  9 ,   B-^  -    -sin  9 ,   Y  =   1  to  an 
expansion  wave  in  which  V  >  1 ,  and  finally  ends  as  a  fast 
wave  beyond  M.  =   1,   B-,  =  sin  0,  V  =  1.   The  limiting  case  0  =  0 
is  not  well  represented  in  this  figure.   It  is  better  plotted 

against  B^  rather  than  r| .   Then  it  appears  as  a  circle  of  unit 

~        2  2 

radius  cent     at  B^  =  0,  M.  =  2,  plus  the  whole  M„-axls.   The 

circle  corresponds  to  the  so-called  switch-on  shocks.   It  is 

possible  for  Bp  to  be  non-zero  in  two  cases;  one  Is  the  0=0 
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case  above  in  which  the  geometry  does  not  distinguish  between 

Bo  and  B-,;  the  other  is  the  Alfven  shock  for  which  M„  -   cos  Q, 
23  A        ' 

~2    ~2      2 
V  =  1,  and  Bo  +  B^  =  sin  9   and  which  marks  the  transition  from 

expansion  wave  to  slow  shocks  in  Figure  2. 

Now  we  do  not  expect  all  these  various  waves  and  shocks 
to  be  physically  significant  in  our  present  problem.   Indeed, 
there  is  some  doubt  as  to  whether  any  of  those  corresponding 
to  points  for  which  B^  <  0  exist  at  all,  since  they  do  not 
satisfy  the  evolutlonarity  conditions  of  Polovin  [29].   Also, 
as  Priedrichs  [10]  shows,  the  expansion  waves  proper  involve 
a  decrease  of  entropy  so  they  cannot  be  significant  when  we 
start  with  an  initially  cold  plasma.   However,  all  the  waves 
do  correspond  to  singular  points  in  the  (B,v)  phase  space  and 
as  such  will  each  affect  the  qualitative  behaviour  of  solutions 
to  some  extent.   Finally,  we  note  that  only  when  M.  <  (7+1)/(y-1^ 
is  there  more  than  one  singular  point  besides  the  initial 
point;  and  this  range  decreases  as  0  increases,  disappearing 
altogether  in  the  transverse  case. 

By  far  the  most  important  singular  points  are  those 
corresponding  to  the  fast  shocks.   For  each  value  of  M.  >  1 
there  is  a  single  such  shock  with  a  "final"  state  for  which 
B^  >  sin  G   and  V  <  1;  and  as  M„  — >  oo,  then  B^  —^   (7+I )  sin  0/(7-1  ^ 
and  V  — >  (7-l)/(7+l).   The  solutions  most  apt  to  the  piston 
problem  should  be  those  joining  the  initial  state  to  this  final 
state,  to  be  identified  as  the  state  at  the  piston.   This  is 
the  case  with  the  one-fluid  theory  where  the  connection  is  by 
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Figure  2.     The  Relationship  Between  M^  and  B3  at  Singular  Points 
Corresponding  to  Shocks,  When  ^=0. 
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means  of  a  simple  jump  discontinuity,  usually  In  all  the 
variables.   But  in  the  two-fluid  theory  this  is  impossible 
since  both  v  and  B  are  always  continuous  functions  of  x 
(except  in  the  limiting  case  where  A  =  O).   The  only  discon- 
tinuities which  are  allowed  are  in  the  fluid  variables,  in 
particular  in  V.   But  V  is  determined  as  a  root  of  the 
quadratic  (3.13)  whose  leading  coefficient  is  a  constant  and 
the  others  continuous  functions  of  B  and  v.   So  V  can  be 
discontinuous  only  by  changing  from  one  root  of  the  quadratic 
to  the  other.   Thus  a  solution  of  the  steady  state  equations 
may  lie  on  one  or  both  of  the  sheets  corresponding  to  the  two 
roots  for  V,  and  a  shock  in  which  V  is  discontinuous  but  B 
and  V  are  continuous  occurs  where  it  jumps  from  one  sheet  to 
the  other. 

At  the  initial  state  the  roots  for  V  are  1  and  (y-1)/(7+1), 
i.e.,  the  initial  state  lies  on  the  sheet  with  the  larger  root. 
Thus  we  shall  call  this  the  "initial"  sheet,  and  the  other  the 
"shock"  sheet.   The  other  singular  points  may  lie  on  either 
sheet.   The  criterion  for  a  state  to  lie  on  the  shock  sheet 
is  easily  seen  to  be 


V  <  -g/2F      or      V  M^  <  ypV  , 


i.e.,  that  the  flow  at  the  state  be  subsonic.   In  terms  of  B^ 

3 

this  takes  the  form 

(3.20)     cos^e  +  (M^-cos^0)B~^sin  6   <  ^7(7+1 ) "^ (2M^+  sin^0-B? ; 
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One  finds  for  the  fast  shocks  that  the  final  state  lies  on 
the  initial  sheet  if  M.  <  Mp(0)  and.  on  the  shock  sheet  if 
M.  >  Mp(0),  where  the  function  n^{9)    varies  from  y2~ at  0=0 
to  2.47  at  9   =  7r/2 .   Its  value  is  marked  with  a  cross  on  the 
curves  of  Figure  2.   The  slow  shocks  also  lie  partly  on  each 
sheet;  at  either  end  of  their  range,  B^  =  +  sin  9,    they  are 
on  the  initial  sheet  but  there  is  a  central  part  of  the  range 
including  the  point  B-,  =  0  which  is  on  the  shock  sheet.   The 
expansion  waves  are  always  on  the  initial  sheet. 

(b )   The  characteristic  roots  at  the  singular  points 
In  linearising  the  equations  to  study  the  solutions  near 
the  singular  points,  the  only  quantity  whose  evaluation  requires 
some  manipulation  is  bv/hB^;    the  most  convenient  form  for  it 
is  found  to  be  -2B^V/(2FV  +  G)  where  F  and  G  are  given  by  (3.l6^ 
Then  the  linearised  equations  for  A  /  0  take  the  form 


'3.21: 


A^=  A  s 
dx 


where  s  denotes  the  vector  (B^ jB^-B^^v^jV^) ,  the  matrix  A  is 
given  by 


(3.22)     A  = 


0    0  0 

0     0  -A 

0  v-iJ,  0 

|i,     0  -p 


and 
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(3.23; 


^i  -  V  -  M^^cos^e  ,      p  =  M'-^Pcos  9 


~2~  A  — 
V  =  2B^V/(2FV  +  G, 


The  corresponding  characteristic  equation  Is 


;3.24)     A^(A^A^  +  p^)  -  A^A(2^L-v)  +  \i{\i-v)    =  0 


2  2       2 
with  discriminant  (Av  +  p  )  -  4Aij.p  .   Thus  there  are  four 

possible  cases  In  the  classification  of  the  roots:  four 

2 
imaginary  roots  from  two  real  negative  roots  for  A  ;  four 

2 
real  roots  from  two  real  positive  roots  for  A  ;  two  real  and 

two  Imaginary  roots;  and  four  complex  roots  when  the  discrim- 
inant is  negative. 

For  the  moment  we  assume  we  are  using  the  mean  scale  length 

for  X  so  that  A  =  1,  P=  e~  -  e.   The  values  of  the  roots  at 

~       ~  _p    2 

the  initial  state,  in  which  V  =  1,  B^  =  sin  G   and  so  v =  M.  sin  9, 

have  been  studied  in  detail  by  Saffman  [33]  who  shows  that  there 

are  a  number  of  distinct  cases.   These  are  indicated  in  Figure  3 

and  are  seen  to  depend  on  whether  f   is  less  than  or  greater 

than  tan  9   and  on  the  value  of  M.  relative  to  M-,  ,  where 

(3.25)    m2  =  1  +  [cos'g(l+r')-l]"  . 

The  U-shaped  curve  given  by  M.  =  M-,  collapses  into  the  e  =  0 
axis  in  the  transverse  case  so  we  always  have  h   real  roots 


-  41  - 


then;  in  the  longitudinal  case,  0=0,  just  the  right-hand 
branch  collapses  into  the  £  =  1  line  so  that  there  are  never 
four  real  roots. 

With  regard  to  the  other  singular  points  a  few  almost 
obvious  facts  are  worth  noting.   Thus,  M-  <  0  for  slow  shocks, 
and  iJ.  >  0  for  fast  shocks  and  expansion  waves;  and  v  >  0  for 
expansion  waves  and  for  fast  and  slow  shocks  when  on  the 
initial  sheet,  while  v<  0  for  shocks  on  the  final  sheet.   The 
change  of  sign  in  v  occurs  where  2FV  +  G  =  0  and  so  entails  v 
becoming  unbounded.   On  the  other  hand,  the  parameter  p  ranges 
over  all  positive  values  independently  of  \i   and  v;  making  it 

sufficiently  small  or  sufficiently  large  will  ensure  that  both 

2 
roots  for  A  are  real. 

For  the  fast  shocks,  the  dependence  of  the  root-type  .on 

e  and  M.  is  shown  in  Figure  4.   One  can  show  (see  Appendix  A) 

that  while  v  >  0,  i.e.,  on  the  initial  sheet,  \x-v   <  0  so  that 

there  are  two  real  and  two  imaginary  roots.   On  the  shock 

sheet  however,  for  each  value  of  M   the  roots  change  in  type 

as  e  increases;  at  e  =  0  they  are  all  imaginary,  at  e  =  1  all 

real  and  for  a  range  of  Intermediate  values  they  are  all 

complex.   This  intermediate  range  collapses  as  9   — >  o  into  a 

o 

single  curve  up  to  M  =  (y+1)/(y-1),  e  =  1;  and  as  0  — >  jr/2 
both  separating  curves  collapse  into  the  e  =  0  axis  so  that  in 

the  transverse  case  all  roots  are  real.   For  slow  shocks  it  is 

2 
clear  from  the  above  observations  that  the  roots  for  A  are 

real  for  all  p .   On  the  initial  sheet  there  are  always  four 
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Figure  3.    The   Nature  of  the  Characteristic   Roots  at  the 
Initial   Singular  Point  for  Q<B<  7r/2,  ^-0, 
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Figure  4.     The  Nature  of  the  Characteristic   Roots  at  the  Fast 
Shock  Singular  Point  for   0<^<7r/2,  P -0, 
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Imaginary  roots,  two  of  which  become  real  at  the  transitions 
to  the  shock  sheet.   But  In  part  of  the  range  on  the  shock 
sheet,  e.g.,  when  B^  =  0,  the  roots  become  all  Imaginary  again. 
Since  the  fast  shocks  and  expansion  waves  are  mutually  exclu- 
sive and  we  need  the  former,  we  shall  not  consider  expansion 
waves  any  further.   This  can  be  accomplished  by  restricting 
our  attention  to  Mach  numbers  greater  than  unity. 


3- 3-   Solutions  In  Special  Cases 

The  theory  of  autonomous  systems  of  ordinary  differential 
equations  Is  particularly  well  developed  for  two-dimensional 
systems;  and  In  that  case  one  can  deduce  a  good  deal  about 
the  solutions  from  little  more  than  a  study  of  the  singular 
points.   It  Is  therefore  valuable  to  study  particular  cases 
of  our  present  equations  In  which  they  reduce  to  only  two 
dimensions.   This  occurs  In  two  distinct  ways  and  we  shall 
study  the  resulting  particular  cases  In  this  section  before 
picking  up  the  general  case  again  In  the  next.   The  termin- 
ology that  we  shall  use  Is  the  same  as  that  of  Nemytskll  and 
Stepanov  [27]  which  can  be  used  as  a  reference  to  the  under- 
lying theory  of  all  of  this  chapter. 

(a)   The  transverse  case 
If  0  =  7r/2  the  magnetic  field  Is  always  transverse  to  the 
direction  of  wave  propagation.   This  Is  the  case  which  has  been 


In  the  general  case  the  lower  limit  Is  (l+p)^,  the  magnetosonlc 
speed. 
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most  extensively  studied  In  the  past  and  thus,  even  though 
it  is  only  a  particular  case  of  (b)  below,  it  deserves 
special  consideration.   However,  the  present  treatment  has 
been  published  previously,  Morton  [26],  so  we  shall  only 
summarise  it  here  for  the  sake  of  completeness  and  as  an 
Introduction  to  the  more  complicated  cases  treater  later. 
More  details  and  references  to  earlier  work  will  be  found  in 
[26]. 

Putting  B-,  =  cos  9  =  0  in  (3-9)-(3-13)  enables  the 
equations  for  Bp  and  v-,  to  be  separated  from  those  for  B^  and 
Vp .   As  Bp  =  v^  =  0  at  all  the  singular  points  we  may  take 
Bp  as  identically  zero  and  consider  only  B^  and  Vp ;  we  shall 
also  drop  the  suffices  in  this  section.   The  equations  are 
completely  independent  of  the  mass  ratio  if  the  mean  scale 
length  is  adopted;*  so  with  A  =  1  they  become 


dB 
dx 


:3.26)  </   g  =  l-BV 


FV^  +  GV  +  H 


V 


There  are  no  slow  shocks  so  that  there  are  always  Just  two 
singular  points.   The  characteristic  roots  of  the  linearised 
system  are  ±{[x-v)^ .      Thus  for  this  system  the  initial  state 


Hence  this  is  a  special  case  of  those  treated  in  the  next 
section  where  we  put  e  =  1  and  leave  0  arbitrary. 
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is  always  a  saddle  point,  and  the  fast  shock  gives  a  centre 

2     2 
on  the  initial  sheet  for  1+p  <  M .  <  M^ ,  and  a  saddle  point 

on  the  shock  sheet  when  M .  >  Mp .   The  same  is  true  for  the 
complete  system  (3.26).   For  the  fact  that  the  equations  are 
analytic  at  the  singular  points  is  sufficient  to  ensure  that 
the  saddle  points  are  retained;  and  the  normally  difficult 
case  of  the  centre  is  simply  resolved  here  because  the  equa- 
tions are  symmetric  across  the  B-axis,  i.e.,  they  are 
unchanged  by  a  simultaneous  change  of  sign  in  x  and  v.   Thus 
the  existence  of  some  periodic  solutions  is  simply  deduced 
for  M  <  M^. 

In  deducing  the  global  picture  we  have  to  take  care  to 
avoid  those  regions  where  V  becomes  complex.   We  describe 
first  the  situation  when  p  =  0  and  take  y  =  2,  as  is  appro- 
priate in  this  transverse  case.   For  M.  less  than  about  1.6, 
V  is  real  for  all  positive  B  (see  Figure  1(a)).   The  solitary 
wave  first  described  by  Adlam  and  Allen  [l]  then  exists  as  a 
loop  in  the  (B,v)  phase  space  about  the  fast  shock  singular 
point.   Inside  it  there  is  a  complete  set  of  periodic  solutions 
of  increasing  entropy  as  the  singular  point  is  approached, 
while  outside,  the  solutions  are  of  negative  entropy  and  thus 
of  no  physical  significance.   This  all  occurs  on  the  initial 
sheet.   The  solutions  on  the  shock  sheet  play  no  role. 

For  somewhat  larger  values  of  M,,,  V  becomes  complex  in  a 
small  region  which  makes  its  appearance  on  the  B-axis  inside 
the  solitary  wave,  as  shown  in  Figure  5(a).   Thus  some  of  the 
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periodic  solutions  no  longer  exist.   When  M.  reaches  2.0  this 
excluded  region  has  swollen  to  such  an  extent  that  it  touches 
the  solitary  wave.   From  this  point  the  solitary  wave  no 
longer  exists,  and  in  the  succeeding  range  2.0  <  M.  <  Mp  -   2.^7 
the  shock  sheet  first  starts  to  be  important.   it  is  then 
possible  to  construct  the  following  solution  which  will  later 
be  seen  to  be  appropriate  to  the  piston  problem.   Starting  from 
the  Initial  point  the  solution  follows  the  incomplete  solitary 
wave  curve  until  by  jumping  to  the  shock  sheet  it  can  follow  a 
trajectory  there  which  touches  the  excluded  region  on  the 
B-axis.   Here  it  can  return  smoothly  to  the  initial  sheet, 
since  the  sheets  meet  on  the  boundary  of  the  excluded  region, 
and  join  the  largest  existing  periodic  solution  about  the  fast 
shock  singular  point.   Such  a  solution  is  shown  in  Figures  5(h) 
and  6.   In  determining  that  such  a  solution  exists  the  chief 
point  that  has  to  be  checked  is  that  the  independent  variable 
decreases  monotonlcally  along  the  whole  solution  curve. 

In  the  final  range  M.  >  Mp  both  singularities  are  saddle 
points  and  are  on  separate  sheets.   There  is  then  a  unique 
trajectory  which  leaves  the  initial  state  jumps  to  the  shock 
sheet  and  enters  the  final  state  on  that  sheet.   In  the  result- 
ing solution,  shown  in  Figure  7^  B  is  a  continuously  differen- 
tiable  and  monotonic  function  of  x,  while  V  has  a  jump 
discontinuity  ahead  of  which  it  is  monotonlcally  increasing 
and  behind  which  it  is  a  monotonlcally  decreasing  function  of  x. 

The  dependence  of  these  transitions  on  p  is  best  shown 
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s\n9** 


Figure  5.    Examples  of  Solutions  inthe  (B,v)  Phase  Plane  for  €*  I 
and  in  the  (Bg.Bj)  Phase  Plane  for  €sO. 

Across  Indicates  a  Saddle  Point,  a  Point  a  Center,  and 
Broken  Lines  Solutions  on  the  Shock  Sheet.     Shaded  Re 
gions  Show   Where    V  is  Complex. 
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Figure  6.     Steady  State  Profiles  for  the  Transverse  Case  Just  After  the 
Solitary  Wave  has  Broken.   9=90°,  M^=2.04,  )9=0.l  . 
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Figure  7.      Steady  State  Profiles  for  the  Transverse  Case  When 

Each  Singular  Point  is  a  Saddle  Point:  ^  =  90**,  M^=2.66, 
^=0.1  . 
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by  plotting  them  against  (l+p)    and  the  compression  ratio  rj, 
which  equals  the  value  of  B  or  l/v  at  the  final  or  fast 
shock  state.   Thus  the  solitary  wave  breaks  when,  at  its 
peak,  the  fluid  speed  equals  the  characteristic  speed;  this 
gives  a  curve  shown  as  Cp  in  Figure  8.   The  other  main  curve 
C^  corresponds  to  M.  =  Mp  and  marks  where  the  second  singu- 
larity becomes  a  saddle  point  on  the  shock  sheet.   Physically 
this  corresponds  to  where  the  flow  in  the  constant  state 
behind  the  shock  becomes  subsonic.   These  two  curves  divide 
the  (t],(1+p)   )  plane  into  three  main  regions:   that  in  which 
the  solitary  wave  exists,  and  those  in  which  solutions  as 
depicted  respectively  in  Figures  6  and  7  exist.   As  p  -^<»the 
two  curves  converge  to  t]  -  l  and  we  approach  ordinary  gas 
dynamics  in  which  all  compression  waves  break  to  form  shocks. 
This  is  the  only  case  in  which  be  whall  consider  the  effect 
of  p .   In  all  the  following  sections  we  shall  assume  p  =  0. 
However,  we  note  that  the  collapse  into  the  one-fluid  model 
as  p.  — >  0  is  quite  general.   For,  if  the  compression  t]  >  i 
then  M.  >  1+p  so  that  for  all  0,e  >  0  the  steady  state 
equations  approach  those  of  the  transverse  case  (apart  from 
a  factor  of  sin  0    in  the  fields)  as  p  — >  oo. 


-  52  - 


(b)   Equal  ion  and  electron  masses 

Putting  e  =  1,  f^  =  0  dissociates  the  equations  for  Bp 
and  B^  as  in  the  transverse  case.   So  we  drop  the  suffices 
on  B-:,  and  Vp  and,  since  the  mean  scale  length  is  again 
appropriate,  put  A  =  1  to  obtain  the  equation 


(3.27)     ^  =  sin  0  -  BV  +  M^^cos^e(B-sln  9\ 


as  the  generalisation  of  the  second  equation  of  (3.26).   The 
characteristic  roots  of  the  linearised  system  are  again  given 
by 

(5.28)    A^  -  (^i-v)  =   0. 

For  0  >  0  the  behaviour  of  the  solutions  is,  except  for 
a  few  modifications,  the  same  as  that  for  the  transverse  case. 
By  utilising  the  simpler  equation  (3.15)  fo^  V  one  can  again 
obtain  an  explicit  equation  for  the  solitary  wave  (see  Saffman 
[35]),  namely. 


(3.29)    B  =  sin  9  +   2m^[sln  0  t  (m^  +  sin^0 ) ^cosh (mx/M^) ]  ^ 


2    2 
where  m  =  M.  -  1  and  the  positive  sign  is  taken.   This  wave 

breaks  when  M.  =  2(1  +  sin  0)    and  the  peak  magnetic  field  is 

2  +  sin  9.      And  when  M.  -   Mp(0)  the  second  singular  point 

switches  to  the  shock  sheet  and  becomes  a  saddle  point.   This 

is  all  a  direct  generalisation  of  the  transverse  case.   For 
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Figure  8.    The   Domains  of  Existence  of  Solutions  to  ttie 
Steady  State  Equtions  in  the  Transverse  Case. 
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7  =  2,  the  differences  start  to  appear  at  0  =  7r/4  and  Tr/6. 
As  we  can  see  from  Figure  2,  when  0  <  7r/4  the  greatest  ratio 
of  the  magnetic  fields  at  the  two  singular  points  Is  no  longer 
3  and  Is  no  longer  attained  as  M.  — >  oo.   instead,  the  field 

ratio  reaches  a  maximum  of  1  +  2  cosec  20  where 

p 
M^  =  (2  +  tan  0)/(l  -  tan  0)  after  which  It  decreases  to  3 

as  M.  — >  00.   The  switch  to  the  shock  sheet  takes  place  before 

this  maximum  Is  reached. 

When  0  <  Tr/6  the  slow  shock  singular  point  Is  present 
even  for  values  of  M  greater  than  unity.   It  Is  a  centre 
when  on  the  Initial  sheet;  and  around  It  there  Is  another 
solitary  wave  given  by  (3.29)  with  the  negative  sign  and 
breaking  when  M.  =  2(1  -  sin  0).   Thus  associated  with  this 
singular  point  Is  a  complete  system  of  solutions  very  similar 
to  those  associated  with  the  fast  shock  singularity  but 
Involving  magnetic  fields  B  <  sin  0.   The  two  systems  are 
well  separated  by  the  Initial  saddle  point  and  can  be  studied 
Independently . 

In  the  limiting  case  0=0  both  shocks  are  called  switch- 
on  shocks,  because  there  Is  no  transverse  field  Initially,  and 
they  differ  only  In  the  sign  of  the  magnetic  field.   In  fact, 
everything  In  the  (B,v)  phase  plane  Is  completely  symmetric 
about  the  v-axls  as  well  as  the  B-axls.   The  relation  between 
Mach  number  and  position  of  the  singular  points  reduces  to 


;5.30)         (M^  -  2)^  +  B^  =  1 
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2 
and  the  switch-on  shocks  exist  only  for  1  <  M.  <  3.   For  higher 

Mach  numbers  the  shocks  degenerate  into  pure  gas  shocks.   We 

note  that  the  three  main  transition  points  in  the  behaviour 

of  solutions  with  increasing  Mach  number  coincide  in  this 

longitudinal  case.   Thus  the  solitary  wave  breaks,  the  shock 

singularity  crosses  to  the  shock  sheet,  and  the  maximum  field 

at  this  singularity  is  attained,  all  at  M,  =  y^. 

(c )   Negligible  electron  mass 

This  case  represents  the  opposite  extreme  to  the  previous 
case  and  together  they  serve  to  bracket  the  general  case.   In 

letting  e  go  to  zero  we  need  to  use  the  ion  scale  length  for  x 

2  2 

so  that  A  =  e  — >0  and  |  '  =  1-e  — >  1.   We  also  have  to  take 

9   <  7r/2.   Then  in  this  case  the  four  dimensional  system  (5-9)- 
(3.12)  reduces  to  two  dimensions  by  (^-H)  and  (3.12)  becoming 
only  algebraic  equations,  i.e.,  by  losing  the  second  deriva- 
tives of  the  fields.   We  get,  after  substitution. 


dR 
;3.31)     M^cos  9.    ^^   =  (M^V  -  cos^0)B^  -  (M^  -  cos^0)sin 


:3.32)     M^cos  0  .  ^=  -(M^V  -  cos^0)B. 


;3.33)     P(B)V^  +  G(B)V  +  H(B)  =  0  , 


where  we  have  explicitly  recognised  that  F,  G,  and  H  now  no 
longer  depend  on  v.   The  characteristic  roots  of  the  system 
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are  given  by  the  equation 


[^.^^)  A^cos^e  +  M^M.([i,-v)  =  0 


which  is  the  correct  limit  of  the  equation  for  the  general 
case  (3.24).   Comparing  this  with  the  equal  mass  case  (3-28) 
we  see  that  for  the  initial  state  and  the  fast  shock  the 
behaviour  of  the  roots  is  completely  reversed.   When  (i-l-v) 
is  positive  there  is  now  a  centre  and  when  negative  a  saddle 
point.   But  for  a  slow  shock,  because  |j.  <  0,  the  type  of  the 
singular  point  is  unchanged.   This  holds  for  both  the 
linearised  and  full  systems,  for  the  equations  still  have 


symmetry  about  the  B^-axis. 


Thus  for  all  M.  >  1  the  initial  point  is  a  centre.   At 
first  sight  this  would  suggest  that  further  analysis  is  point- 
less since  no  solution  can  then  arise  from  the  initial  state. 
However,  various  modifications  to  the  equations,  e.g.,  the 
introduction  of  slow  time  dependence,  can  allow  the  initial 
point  to  be  approached  and  then  the  present  solutions  are 
important  approximations.   The  fast  shock  point,  on  the  other 
hand,  is  a  saddle  point  on  the  initial  sheet  for  1  <  M.  <  Mp(0) 
and  becomes  a  second  centre  when  it  shifts  to  the  shock  sheet 
at  Mp(0).   So  for  sufficiently  small  M.  we  expect  there  to  be 
a  solitary  wave  beginning  and  ending  at  the  fast  shock  point. 
Except  when  0=0,  this  is  more  difficult  to  demonstrate 
explicitly  than  in  the  equal  mass  case  because  the  full  quadratic 
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for  V  has  to  be  used.   However,  we  can  prove  Its  existence 
for  certain  values  of  M,  as  follows: 

For  some  a  >  0  to  be  chosen  later,  let 


R^  =  Bg  +  (B^  -  a  sin  9)    ,      Bg  =  R  cos 


Then  from  (3-31)  and  (3-32)  ee  obtain 


2 
i  COS  9  •  ^^  =    -  M  [1  -  aV  -  M^  cos  0(l-a)]B2Sin 


(3.35) 


i.e.   f^  =  -  K  cos 
dx 


For  0  <  0  <  ir/2 ,   we  have  K  >  0  if  a  is  chosen  so  that 


V<a^-(a^-  l)M^^cos^( 


But  on  the  trajectory  leaving  the  fast  shock  point  the 
pressure  is  always  positive;  thus  from  {'^•^)    and  (3.5) 

M^(l-V)  >  -  I  sin^e 


]   -P   1   -P    2     "5    1-2    2 
I.e.  V  <  1  +  ^  M^^-  i  M^  cos  0  <  ^  -  ^  M^  cos  ( 


so  that  the  required  inequality  is  certainly  satisfied  for 
a  =  2/3. 

We  can  now  prove  the  existence  of  the  solitary  wave  in 
those  cases  where  (i)  the  fast  shock  point  is  on  the  initial 
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sheet,  and  (11)  the  value  of  V  on  the  Initial  sheet  Is  real 
and  positive  everywhere  Inside  the  circle  C  centred  at 
Bp  =  0,  B^  =  a  sin  9   and  passing  through  the  fast  shock 
point.   We  know  that  there  are  such  cases.   For  when  0   and 
M.  are  sufficiently  small  V  Is  never  complex;  also,  on  the 
Initial  sheet.  Re  V  >  0  In  the  disc  G  <  0,  which  can  easily 
be  shown  to  cover  the  circle  C  for  all  a.   In  general,  we 
have  shown   that  V  Is  complex  In  a  roughly  annular  region  about 
a  central  core,  In  which  V  >  0  and  which  contains  the  whole  of 
the  B^-axls  from  the  origin  to  the  fast  shock  point. 

So  now  consider  that  part  of  the  right  half-plane  bounded 
by  the  circle  C,  any  one  of  the  periodic  solutions  about  the 
Initial  point,  and  the  B^-axls  (see  Figure  9)-   It  Is  clear 
that  the  trajectory  from  the  fast  shock  point  enters  this 
region.   But  there  are  no  singular  points  or  limit  points  In 
the  Interior  because  dR/dx  <  0  there;  for  the  same  reason  the 
trajectory  cannot  leave  by  crossing  C.  Hence  It  must  cross  the 
B^-axls,  and  closing  It  up  by  symmetry  yields  the  solitary 
wave.'   This  must  Include  a  singular  point  and  hence  encircles 
the  Initial  state  point. 

The  solitary  wave  will  first  cease  to  exist  for  some 
Mach  number  not  greater  than  Mp(0).   The  solutions  In  the 


See  9  3»1  and  appendix  A. 

'   The  possibility  that  the  trajectory  enters  a  slow  shock 
singularity  can  be  discounted  on  grounds  of  entropy 
difference . 
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intervening  cases  up  to  M„(0)  will  probably  behave  as  in  the 
similar  intermediate  cases  of  section  (a)  above;  but  to 
ascertain  this  would  require  a  numerical  investigation  which 
we  have  not  carried  out.   When  M»  >  M^{0)    the  periodic  solu- 
tions about  the  fast  shock  point  are  limited  by  the  annulus 
of  complex  V-values.   This  defines  a  largest   one  which  will 
touch  the  annulus  and  will  also  touch  a  periodic  solution 
around  the  initial  point.   Together  these  define  a  complete 
solution  which  has  a  jump  discontinuity  in  V  where  they 
touch  (see  Figure  5c). 

In  the  limiting  case  0=0  matters  are  somewhat  simpler. 
It  follows  from  (5.35)  that  all  trajectories  are  circles  about 

the  initial  point,  the  origin;  they  are  pure  Alfven  waves 

2 
with  V  constant  on  each.   For  M  <  2  they  lie  inside  a  circle 

of  radius  given  by  (3.^0)  every  point  of  which  is  singular. 

These  points  are  the  limits  of  the  slow  and  fast  shock  points. 

Outside  this  circle  lies  the  annulus  in  which  V  is  complex; 

2 
it  approaches  as  M.  increases  to  the  value  2,   when  its  inner 

boundary  is  coincident  with  the  circle  of  singular  points. 
As  M.  increases  further,  it  recedes  and  is  filled  by  the 
circular  trajectories,  for  the  singular  points  are  now  on  the 
shock  sheet.   But  in  this  case  we  do  have  complete  solutions 
which  look  a  little  nearer  those  required  for  the  piston  prob- 
lem.  Each  singular  point  on  the  shock  sheet  is  a  constant 
state  to  be  identified  with  that  at  the  piston;  and  it  can  be 
directly  connected  by  a  Jump  to  the  initial  sheet,  accompanied 
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Figure  9.    Diagram  for  Proving  the  Existence  of  a  Solitary  Wave 
in  the  Case  of  Negligible  Electron  Mass. 
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by  a  discontinuity  In  V,  to  a  sine  wave  oscillation  about 

2 
the  Initial  state.   Note  that  the  critical  value  M,  =  2  refers 

to  the  case  y  =  2. 

3.^.   Solutions  In  the  General  Case  .:..;.o 

It  Is  valuable  to  begin  by  briefly  reviewing  the 
information  collected  in  the  two  previous  sections  on  the 
singular  points  and  the  special  cases.   The  equations  (3.9)- 
(3- 13)  describing  the  steady  flows  depend  upon  three  parame- 
ters, e,    Q   and  M.,  but  two  of  their  most  important  features 
do  not  depend  upon  the  first.   Thus  if  the  mean  scale  length 
is  used,  the  quadratic  determining  the  specific  volume  is 
Independent  of  e.   The  same  is  true  of  the  position  of  the 
singular  points,  including  the  value  of  M.  at  which  the  fast 
shock  singularity  shifts  to  the  shock  sheet.   Hence  the 
overall  geometry,  of  the  (Bj,v)  phase  space  depends  only  on  0 
and  M. .   In  the  special  case  e  =  1  the  phase  plane  given  by 
Bp  =  V-,  =  0  was  studied,  while  in  the  case  e  =  0  that  given 
by  Vp  =  v^  =  0  had  to  be  used.   These  two  extreme  cases,  in 
fact,  serve  to  cover  this  facet  of  the  problem  very  well. 
The  type  of  the  singular  points,  however,  and  therefore  the 
nature  of  the  solutions,  depend  markedly  on  e  and  are  very 
different  at  the  two  extremes.   If  we  consider  only  the  profile 
of  B^,  it  always  has  the  general  form  of  a  smooth  wave  front. 
But  for  the  equal  mass  case  the  front  always  rises  exponentially 
from  the  initial  state  and,  when  M.  <  Mp(0),  is  succeeded  by 
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an  oscillatory  wave  train,  while  when  e  =  0  there  Is  always 
an  oscillatory  precursor  to  the  front  and  oscillations  behind 
only  when  M.  >  n^{9) . 

The  main  object  of  study  In  this  section  is  how  this 
changeover  is  brought  about  as  e  varies  from  unity  to  zero. 
We  shall  not  give  as  many  quantitative  results  as  in  the 
special  cases,  although  these  could  easily  be  obtained  by 
numerical  integration  of  the  equations,  for  it  is  sufficient 
to  describe  the  qualitative  behaviour  at  the  singular  points 
and  its  dependence  on  the  parameters  e,    0   and  M..   This  will 
be  done  only  for  the  linearised  equations  (5-21).   Most  results 
extend  to  the  complete  equations,  the  only  doubtful  phenomenon 
being  the  generalised  vortex  which,  however,  is  less  important 
than  the  centres  occurring  in  the  two-dimensional  systems. 

The  eigenvalues  A  of  the  matrix  A  of  (5.22)  are  given  by 

(3.36)  2A^A^  =  A(2n-v)  -  p^  +  [  (AV  +  p^f    -    4A^P^]2  . 
The  corresponding  eigenvectors  are 

(3.37)  B2  :  By   B^  :  V2  :  v^  =  pA  :  AA^-^l  :  -A(AA  -^i)  :  pA   . 

Let  us  consider  first  the  case  where  e  is  near  unity;  then  we 
take  A  =  1,  and  p  «  1  to  get 

(3.38)  a2  =:  n(l  -  p2/v)  ,   (H-V)(l  +  p^v)  . 

There  are  two  cases,  M.  <  M^{9)    and  M  >  M^{9): 

(i)   e  ^  1,  M.  <  M2(9) •   The  initial  point  is  a  generalised 
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saddle  point  with  four  real  eigenvalues,  two  positive  and  two 
negative.   The  eigenvectors  corresponding  to  the  latter  yield 
a  two -parameter  family  of  solutions  leaving  the  initial  point. 
With  these  projected  onto  the  (Bp,B-,)  plane,  the  situation 
looks  very  similar  to  that  for  a  node  in  a  two-dimensional 
system;  all  the  solutions  leave  the  initial  point  along  the 

eigenvector  corresponding  to  the  eigenvalue  which  is  smallest 

-  2    - 

in  absolute  value,  i.e.,  {\i-v)^{l  +   p  /v)^.   From  (3.37)^ 

this  implies  that  at  the  front  of  the  progressing  wave  the 
magnetic  field  rises  exponentially  from  its  initial  value  with 

(3.39)  Bp  :  B^-sin  9   =    (e"^-  e)(M^-  l)^cos  0  cosec^S 

where  we  have  substituted  for  iJ.,  v,  and  p  from  (3-23).   This 
formula  shows  how  the  Bp  component  increases  as  we  leave  the 
equal  mass  case.   And  it  indicates  the  complexity  at  0  =  0, 
e  =  1,  the  limits  along  8=1  and  along  0=0  being  entirely 
different  in  nature. 

At  the  fast  shock  singularity,  we  have  v  >  |j.  >  0  so  that 
there  is  a  composite  saddle  point  with  two  imaginary  eigen- 
values and  two  real,  of  which  one  is  negative  and  one  positive, 
Corresponding  to  the  last  is  an  eigenvector  with 

(3.40)  B^  :  (B^  -  B^)  =  -V  :  pH* 

while  the  imaginary  roots  give  a  periodic  solution  for  which 

(3.41)  Bp  :  (B^  -  B^)  =  tip(v-n)*  :  V  . 
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Thus  behind  the  wave  front  It  is  mainly  the  B^  component  that 

Is  oscillatory  while  Bp  Is  exponentially  damped. 

(11)  e  r^  1 ,   M„  >  Mp(0).   The  situation  at  the  Initial  point 

Is  unchanged  and  at  the  fast  shock  point,  since  v  <  0  <  |i.-v, 

the  situation  Is  now  similar.   The  final  state  Is  approached 

exponentially  with  the  limiting  field  direction  given  by 

(3.40). 

At  the  other  extreme,  where  e  is  small,  we  need  to  employ 

2 
the  ion  scale  length  so  that  p  remains  bounded.   Then  A  =  e  , 

I   =1-6   and  we  get 

2A^A^  =  A(2^L-v)  -  p^  +  [p^-  A(2^L-v)  -  ^   A^p~^(2^-v)^ 

(3.42)  +  I  A^p-^v^  -...],      .... 

2       2-4  -2 

I.e.   A   =  -  p  £   ,     -|i,(|a-v)p 

In  all  cases  the  first  pair  yields  oscillations  whose  period 

o 

tends  to  zero  like  e   and  which,  as  we  see  from  (3.37)?  is 
symmetric  in  Bp  and  B-, .   This  adds  a  high  frequency  component 
to  all  solutions.   We  also  note  that  if  M.  =  0(e   )  then  the 
period  has  the  scale  of  the  mean  ion-electron  gyro-radius. 
Saffman  [32]  has  shown  that  a  solitary  wave  with  these 
characteristics  exists  for  all  values  of  e  when  0=0. 

For  the  other  pair  of  roots  we  have  the  usual  two  cases: 
(ill)   e  .-^  0,  M»  <  Mp(0)  .   At  the  initial  point  these  roots 
are  also  Imaginary,  giving  oscillations  with  a  period  of 
approximately 
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(5.43)  27ri^^[ (M^  -  1)(M^  -  cos^0)]~2cos  9 
and  elliptic  polarisation  given  by 

(3.44)  B2  :  B^  -  sin  9   =  ±    1 (M^  -  l)^  :  (m^  -  cos^e)^  . 

At  the  fast  shock  singularity  the  roots  are  real  and,  apart 
from  the  high  frequency  oscillations,  there  is  exponential 
convergence  to  the  final  state  along  Bp  :  (B^-B^)  =  (v-|-l)^  :  [x' 
(iv)   e  ~  0,  M.  >  Mp(9).   Here  the  only  change  is  that  the 


roots  at  the  fast  shock  singularity  are  also  imaginary  so  that 
there  are  oscillations  before  and  behind  the  main  wave  front. 

A  glance  at  Figures  3  and  4  shows  that  when  M. '~-  1  the 
situations  described  under  (i)  and  (iii)  above  cover  the  whole 
range  of  e;    the  dividing  line  occurs  where 

(3.45)     e"^  -  £  =  tan  9 

which  agrees  with  the  deductions  drawn  from  the  asymptotic 
theory  of  j2.3-   At  higher  Mach  numbers,  regions  where  the 
eigenvalues  are  complex  intervene.   These  give  damped  oscilla- 
tions and  the  behaviour  of  solutions  is  always  more  similar 
to  the  e  -^  1  cases  than  for  the  lower  Mach  numbers. 
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3.5-   The  Effect  of  Friction 

Finally,  in  concluding  our  study  of  steady  flows  we  show 
that  the  introduction  of  friction  will  sometimes  allow  the 
initial  and  fast  shock  singularities  to  be  joined  by  a  single 
continuous  solution.   The  postulation  of  friction  forces 
between  the  motions  of  ions  and  electrons  results  In  the 
Introduction  of  a  term  cc J  on  the  right-hand  side  of  equation 
(2.5)  (see,  e.g.,  Morawetz  [24]).   This  causes  terms  kv  to  be 
added  to  the  right-hand  sides  of  equations  (3-11)  and  (5-12), 
where  k  >  0;  and  in  fact   this  is  the  only  change  to  the 
steady  state  equations. 

We  note  first  that  this  turns  all  the  centres  and 
generalised  vortices  into  foci.   For  the  characteristic  equa- 
tion for  the  linearised  system  becomes 


0A6)  (AA^  -  kA  -  ^L)(AA^  _  kA  -  ^i  +  v)  +  p^A^  =  0 


In  the  equal  mass  case,  p  =  0  and  every  root  is  prevented  from 
being  pure  imaginary  by  a  positive  real  part  ^   k/A.   Thus  the 
oscillations  behind  the  wave  front  are  damped.   At  the  other 
extreme  the  equation  becomes 

(3.47)     (p^  +  K^)A^  +  k(2m.-v)A  +  n(n-v)  =  0  . 

To  obtain  a  non-zero  imaginary  part  to  an  eigenvalue  requires 
l-L-v  >  0  so  that  K  contributes  a  negative  real  part.   This  again 
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causes  damping  of  the  precursor  away  from  the  wave  front.  To 
check  Intermediate  cases  we  substitute  a  pure  Imaginary  root, 
la.  Into  {'^A6)    to  obtain  the  condition 

K(2ia  -  V  +  2Aa^)  =  0  . 

Then  substitution  of  this  expression  for  v  Into  {^A6)    yields 

(Aa^  +  \if  +  o^{k^   +  p2)  =  0 


which  can  clearly  never  be  satisfied. 

Now  to  prove  the  existence  of  solutions  joining  the 

singularities,  we  consider  first  the  behaviour  of  the  entropy 

*  2 

variable  A  =  pV   on  the  Initial  sheet.   Substituting  for  p 

yields  various  expressions  for  A: 
(3.^8a)    A  =  -V^(^  G  +  M^V) 

(3.^8b)    A  =  -(^)[G  +  2(G^-  4fh)*] 


OASc)        A   =   |v^[sln^0  -  Bg  -  b|  +  2M^(1-V)]  . 

Thus   A   Is  a  continuous  real  function  of  the  variables  Bp , 

2       2     2 
B^,  1^,  where  as  before  ^   =  A(vp  +  v^),  as  long  as  V  Is  real. 

Prom  (3-'4-8b)  It  Is  clear  that  the  condition  A  >  0  Implies 

G  <  0  and  3G  <  16fH.   The  first  defines  the  cylinder  H,    and 

the  second  a  torus  J    which  Includes  the  region  where  V  Is 

complex  (see  Figure  9)  and  whose  central  "hole"  consists  of 

a  single  point  --  the  Initial  singularity  where  A  ==  0. 


We  assume  y  =  2  In  this  section. 


Further,  It  turns  out  that  the  conditions  for  a  stationary 
value  of  A   In  this  3-dlmenslonal  phase  space  are  just  those 
for  a  singular  point.  I.e.,  (;5.  l6)  -  (3  •  l8)  .   For  A  depends 
quadratlcally  on  Bp ,  Vp ,  and  v^  so  these  must  all  be  zero  at 
such  a  point,  while  for  B^  we  have 

ir=  -v(|g+3m2^)(||-)  -b/ 

=  -V[sln  e  +   M~^cos^0(B^-  sin  0)  -  2VB^]  -  B^V^ 

=   M~^V[ (M^V  -  cos^0)B^  -  (M^  -  cos^e)sln  0]  . 

Finally,  the  energy  conservation  equation  (5-1)  Is  unchanged 
by  the  Introduction  of  friction  and  yields  for  A  the  equation 

a^A  =  kV^^  . 

Along  any  trajectory  of  the  steady  state  equations  this  becomes 

(3-49)    ^  =  -  kM-^^^  . 

Let  us  now  restrict  ourselves  to  those  cases,  considered 
previously  In  p3-5(c)  for  which  V  Is  complex  at  no  point  of 
the  phase  space.   Then  ^  ^  J    so  that  V  Is  real  and  positive 
everywhere  on  the  Initial  sheet  Inside  J  .      As  we  have  seen 
above  there  Is  a  trajectory  leaving  the  Initial  point  with 
decreasing  x  and  thus,  by  (3-^9)^  with  Increasing  A.   It 
must  therefore  enter  J    from  which  It  cannot  escape,  since 
A  =  0  on  the  boundary.   Inside  j  It  can  be  continued 


indefinitely  until  it  enters  a  singular  point,  any  other 
result  being  precluded  by  (3-^9) •   The  only  remaining  question 
is  whether  it  enters  the  fast  or  slow  shock  singularity.   In 
the  particular  case  of  equal  ion  and  electron  masses  this  is 
simply  answered.   For  then  solutions  are  restricted  to  the 
section  of  the  torus  ^7  given  by  setting  Bp  =  0.   The  boundary 
curve  consists  of  the  two  solitary  wave  lobes,  one  about  the 
fast  and  one  about  the  slow  shock  singularity,  which  have  only 
the  initial  point  in  common.   Thus  a  trajectory  entering  one 
cannot  stray  to  the  other  and  there  are  solutions  joining  the 
initial  point  to  each  of  the  shock  points.   It  is  also  clear 
that  in  this  case  one  can  extend  the  present  result  to  larger 
values  of  M„  and  9,   where  V  may  be  complex  at  some  points  of 
J      so  long  as  these  do  not  lie  in  the  particular  solitary 
wave  lobe  being  considered.   For  the  fast  shock  point  this 
covers  a  non-empty  range  of  M.  for  every  value  of  9,    ranging 

from  l<M^<1.5at    0   =   0  to   1<M^  ^2. 5  at   0   =  tt/2. 

2 
For  general  values  of  the  mass  ratio  e  ,  it  is  much  more 

difficult  to  demonstrate  that  there  is  a  trajectory  leaving 

the  initial  point  and  entering  the  fast  rather  than  the  slow 

shock  point.   For  when  V  is  real  everywhere  in  J  ,  the  slow 

shock  point  must  be  present  and  the  absolute  maximum  of  A 

must  in  fact  be  attained  there.   Otherwise  the  maximum  would 

2 
have  to  be  attained  at  the  fast  shock  point.   But  for  some  e 

there  is  a  trajectory  leaving  this  point  with  constant   A,  and 

points  along  this  trajectory  are  not  stationary  points  except 
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2 
In  the  special  case  0=0  when  M„V  =  1  there.   Thus  In  general 

the  fast  shock  point  Is  a  saddle  point  for  A  and  the  maximum 

occurs  either  at  the  slow  shock  point*  or  on  the  surface  where 

V  becomes  complex. 

2 
Finally,  we  note  that  at  the  other  extreme,  e   =  0,  this 

question  Is  again  easily  resolved  when  the  conditions  hold 

under  which  the  existence  of  the  solitary  wave  was  shown  In 

§  3-3(c).   For  then  solutions  lie  on  the  plane  ;;  =  0  and  there 

Is  a  trajectory  leaving  the  fast  shock  point  and  entering  the 

solitary  wave  loop  with  decreasing  A.   By  arguments  Identical 

with  those  above,  It  Is  easily  seen  that  Is  must  enter  the 

Initial  point. 


The  slow  shocks  occurring  here  are  called  "Improper"  by  Bazer 
and  Erlcson  [5]  who  showed  that  their  entropy  Increase  was  of 
the  first  order  In  the  shock  strength  rather  than  the  usual 
third  order. 
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4 .   Solution  of  the  Piston  Problem 

^ . 1 .   Description  of  the  Problem 

We  shall  here  consider  the  solution  by  finite  difference 
methods  of  the  following  problem.   A  cold  plasma  Initially 
occupies  the  half  space,  x  >  0,  bounded  by  a  fixed,  plane, 
conducting  wall  at  x  =  0.   The  plasma  is  at  rest  and  contains 
a  uniform  magnetic  field  of  unit  intensity  and  inclined  at  an 
angle  9   to  the  x-axls.   Either  an  electric  or  a  magnetic 
field  is  now  applied  at  the  wall  and  in  its  plane,  and  raised 
to  a  constant  final  value.   As  a  result,  the  plasma  separates 
from  the  wall  leaving  behind  it  a  vacuum  region.   The  motion 
of  the  plasma  is  assumed  to  satisfy  the  equations  of  |2.1 
with  the  origin  of  the  Lagranglan  variable  x  taken  at  the 
plasma -vacuum  interface. 

The  Initial  conditions  are  clear.   For  the  boundary 
conditions  at  x  =  0,  the  situation  when  9  f  it/2    is  rather 
different  from  the  special  case  9   -   Tr/2 .   As  pointed  out  by 
Prledrlchs  [10],  only  in  the  latter  case  can  a  contact  discon- 
tinuity be  maintained.   Thus,  in  general,  the  single  condition 
on  the  fluid  variables  which  the  characteristics  Indicate  is 
needed,  will  be  density  p(0,t)  =  0.   When  9   =   ■k/2   this  could 
be  replaced  by  the  condition  of  total  pressure  balance,  but 
for  uniformity  we  shall  continue  to  use  p(0,t)  =  0.   Under 
this  condition  the  tangential  electric  and  magnetic  fields 
will  be  continuous  across  the  Interface.   Since  there  is  no 
current  in  the  vacuum  the  magnetic  field  is  uniform  there  and 
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we  denote  it  by  B^'^'^(t).   Then  when  the  Interface  is  distance 
X  from  the  wall  the  electric  field  there  will  be 
e(^)  +  X(i  X   ^Js^^h   where  E^^^  is  the  field  at  the  wall  and 
i  Is  a  limit  vector  in  the  x-direction.   Thus  we  have  two 
alternative  boundary  conditions  depending  on  whether  the 
magnetic  or  electric  field  is  given  just  inside  the  vacuum 
at  the  wall: 

(4.1)  p(0,t)  =  0  J 
and 

(4.2)  B(0,t)  =  B^''^(t) 
or 

(4.5)      E^^^  +  Xi  X  S^B  =  B  y  u  +  Al  ><  ^t^x^  ^  r^B-L^^B  , 

where  equations  (2.l6)  and  (2.17)  have  been  used  to  obtain 

(4.3)  and  the  arguments  suppressed  for  the  sake  of  clarity. 
The  distance  X  is  given  by  integrating  dx/dt  =  u^(0,t). 
After  integration  with  respect  to  t  the  last  equation  takes 
the  more  useful  form 

(4.4)  Al  X  S^B  -XixB+B^    (ixu  +  r^x^)dt  = 


E^^^t)dt 


Note  that  only  the  y-  and  z -components  of  (4.2)  and  (4.4) 
are  meant  to  be  used,  giving  Just  two  boundary  conditions  on 
the  fields.   This  is  the  correct  number  in  all  cases;  for, 
although  in  the  limiting  case,  A  =  0,  the  order  of  the  equa- 
tions (2.14),  (2.15)  Is  reduced,  the  characteristics  which 
which  are  lost  are  parallel  to  the  t-axis. 
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H-  .2  .      Finite  Difference  Equations 

Second  order  accuracy  Is  maintained  In  the  difference 
equations  by  using  a  "staggered"  mesh.   With  a  constant  mesh 
size  h  and  a  time  step  k,  u  Is  correctly  centred  at  points 
X  =  jh,  t  =  (n-l/2)k,  (j,n  =  0,1,...),  where  Its  value  Is 
denoted  by  u.  with  components  u.  . ,  1  =^  1,2,3.   Similarly 

V,  p,  p.  A,  B  are  all  centred  at  x  =  (j  +  l/2)h,  t  =  nk, 

n 
where  their  values  are  denoted  by  V.  etc.   Finally,  to  allow 

Integration  of  the  equations  through  shocks,  an  artificial 

viscous  pressure  q  Is  Introduced  In  the  manner  suggested  by 

von  Neumannand  Rlchtmyer  [28];  this  Is  centred  like  p  but  the 

product  qV  Is  centred  at  x  =  (j  +  1/2 )h,  t  =  (n  -  l/2)k  and 

again  denoted  by  (qV)..   With  this  notation  and  writing 

D  [w.l  for  h~  (w . , -,  -  w.),  the  equations  become 
+   J         ^  J+1    J^ 

n,    n\  n+1      n       ,^  r  n  ,   n  ,  l/^^n   >,2    l/'-^n  \2 

(^.5)      u^^  .^^  =  u^^^  .^^  -  kD^[p.  4-  q.  +  2(B2^  .)   +  ^{B^ ^  .) 


with  p  =  A/y'^   =  Ap'^ 


(4.6)      V"+^  -  V"}  +  kD^[u'?+^]  , 


(^•7)  (qV)5+^  =  a^h^  ^mln(0,D^[u^+J])]^ 

•     •  n+1        1,  /   ,,  sn+l/^^n+l    ,    -t^ns-1  n 

giving  qj        =   4(q/)^.       (V^.        +  ^ j)         "   ^j    ' 

/I,    ON  n^+l         .n         /      -,  w   ^^Nn+1/1  ,,n+l    ,     1  ,,nsY-2,  _    r    n+1-, 

(4.8)  A  =A      -    (Y-l)(qV).       (2  V  +2^1^         kD^[u^     .]     , 
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and 
(^•10) 


^rVi  =  ^?  -14-1  +  kBiD^[By  A    ,         1  -  2,3  , 


-a  .  B  . ,  -,  +  b  .   B  .    -  c  .  B  .  -,  =  d  . 


In  the  last  equation 


-  I  kfE^    A 


and 


2  n 
2A  +  h  V  . 

J 


-kPB^ 


2..n 


kTBi 


2  n 

2A+  h  V. 

J 


/ 


where 


d"?  =  (hn^-?-A5^  )B^  +  kh^B,  D^ 
J    ^    J       .J       1  + 


6w.=:w.-,  -2W.  +  W.,-, 
J     J-1      J     J+1 


n+1' 

n+1 
5, J 


IkTB^ 


-6 


B^ 
J 


The  values  of  all  the  variables  at  n  =  0  are  prescribed 
by  the  initial  conditions  and  passage  from  time  step  n  to 
(n+l)  is  accomplished  by  solving  equations  (4.5)  to  (4.10)  in 
the  order  given.   The  index  j  normally  runs  from  zero  to  some 
integer  N(n).   There  is  some  complication  at  the  plasma- 
vacuum  interface  caused  by  the  boundary  conditions,  principally 
the  condition  p  =  0  or  V  =  <».   However,  because  of  the  mesh 
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adopted  the  definition  of  V  at  the  Interface  Is  avoided,  so 
(4.6)  can  be  used  even  for  j  =  0  although  It  Is  preferable 
to  use  the  corresponding  equation  for  p  when  j  Is  small.   Also 
B.  Is  determined  by  (4.10)  back  to  j  =  -1,  so  that  u   Is  given 
t>y  (4.5)  and  (4.9)  by  setting  D  [p  +  q  ]  =  0.   To  complete 
the  system  (4.10)  the  boundary  conditions  on  the  fields  are 
used.   If  the  magnetic  field  Is  given  then  we  take 


:4.11)     B^^  +  ^S  =  2B^''^(nk; 


and  If  the  electric  field  Is  given  (4.4)  Is  used  with  B 
replaced  by  i(B^.  +  b"),   S  B  by  h~^(B"  -  B^J  and  the 

(—      J-        O         yC  O        ~_L 

Integrals  evaluated  by  the  trapezoidal  rule.  The  resulting 
system  of  equations  Is  block  trl-dlagonal  and  may  be  solved 
by  the  usual  recurrence  relations  as  follows:   E  and  P  are 

'^  0         0 

determined  by  comparing  the  boundary  conditions  with  the  equa- 
tion B  -I  =  E  B  +  F  ;  then  after  a  little  simplification  we 
-1    o  o    o'  ^ 

find 

E.^,  =  (aT^b  .  -  E.)"^ 

(4.11)  J  =  0,1,2,. ..,N    ^    ■ 


-1^ 

j+1    "j+l^"j  ^j  '  "j 


P-. ,  -,  =  E,,_^(a  .  d  .  +  F 


followed  by 


(4.12)     B^^,  =  B(t=0),   B.  =  E.^^B.^^+F.^^  , 


N,N-1, . . . ,0,-1  , 
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where  the  superscripts  have  been  suppressed.   The  value  of  N 
Is  determined  at  each  time  step  by  continuing  (^.11)  until  It 
Is  assured  that  B,,  Is  sufficiently  close  to  B^_^-^-      Because 
the  mean  scale  length  (in  which  A  =  l)  Is  always  used  when 
9   =  7r/2,  the  matrix  a.  Is  always  non-singular.   For  the  same 
reason  we  can  deduce  from  the  observation  I|e  ||  <  1^  that 
I|E.  II  <  1  for  all  j  >  0.   Thus  the  above  process  Is  numerically 
very  stable  In  both  the  forward  and  the  backward  recursion. 

The  stability  and  convergence  properties  of  the  finite 
difference  scheme  as  a  whole  are  discussed  In  Appendix  B. 
When  no  artificial  viscosity  Is  needed.  I.e.,  when  no  breaking 
occurs.  It  Is  shown  that  a  sufficient  condition  for  the  solution 
to  converge  to  the  correct  solution  as  h  — >  0  Is  that 

(4.15)     (k/h)c    <  1 
^    -^  ^  ^   I     '    max 

*  ,    2  2  \-k 

where  c  Is  the  magneto-sonic  speed  (c„  +  Ca)  •   This  result 

Is  obtained  by  methods  based  on  the  work  of  Krelss  [l8]  using 

energy  Identities.   In  the  present  case  the  differential 

equations  yield  the  following  energy  identity: 

(4.14)     [Pu  +  (E  +  u  /  B)x  B^]  .  1 


I 

c» 

r 


_d_ 
dt 


[^u-u  +  #(S  B-S  B)  +  (y-1)"^pV+ J  VB^-B'^ldx 

£-         ^   X    X  ^ 


o 


where  the  left  hand  side  is  evaluated  at  the  plasma-vacuum 
interface  and  P  denotes  the  total  pressure.   Besides  providing 


The  speeds  here  relate  to  the  Lagrangian  co-ordinate;  see  the 
appendix. 
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the  motivation  for  the  energy  identities  for  the  difference 
equations,  this  identity  provides  a  valuable  check  on  the 
accuracy  of  the  computations. 

When  shocks  form  and  the  artificial  viscosity  is  needed 
there  is  a  further  restriction  on  the  time  step.   This  is 
determined  by  the  behaviour  in  the  narrow  shock  region,  where 
the  fields  are  almost  constant  so  that  the  restriction  is  the 
same  as  that  for  ordinary  fluid  flow.   A  stability  condition 
for  this  case  is  given  by  Richtmyer  [30],  which  for  7  =  2  is 

where  r\   is  the  compression  ratio  of  the  shock.   This  was 
found  to  be  quite  satisfactory  in  practice. 


■^.3-   Summary  of  Results 

As  a  result  of  the  considerable  amount  of  qualitative 
information  obtained  about  the  steady  flows  it  is  necessary 
to  compute  only  relatively  few  time  dependent  cases.   Such 
computations  at  a  few  selected  points  in  the  parameter 
space  give  qualtltative  results  which  are  in  such  good  agree- 
ment with  expectation  that  the  whole  picture  can  be  readily 
filled  in. 

Although  solutions  depend  to  some  extent  on  whether  an 
electric  or  magnetic  field  is  imposed  at  the  wall,  and  on  its 


orientation,  these  effects  are  mainly  confined  to  the  neigh- 
borhood of  the  plasma -vacuum  interface.   Thus  the  main  body 
of  results  have  been  obtained  by  imposing  a  magnetic  field  in 
the  z-direction.   For  an  initially  cold  plasma,  then,  these 
results  depend  on  three  parameters  as  shown  in  Figure  10:  the 
initial  field  orientation  9,    the  mass  ratio  e  ,  and  the  Alfven 
Mach  number  M,  (or  the  final  strength  of  the  imposed  field). 

In  the  two  limits  (i)  9   =   7r/2 ,  0  <  e  <  1,  and  (ii)  e  =  1, 
0  <  0  <  ir/2 ,    a  steady  flow  satisfying  the  initial  conditions 
exists  when  M„  >  2(1  +  sin  0),  i.e.  after  breaking  occurs. 
Then  the  time  dependent  solution  rapidly  converges  to  this  as 
was  shown  in  the  earlier  paper  [26]  on  the  transverse  case. 
Even  for  smaller  Mach  numbers  the  time  dependent  solutions  can 
be  quite  accurately  predicted;  the  amplitude  of  the  wave  front 
and  the  oscillations  behind  it  can  be  obtained  from  the  soli- 
tary wave  solution  (3-29);  and  the  period  of  the  oscillations 
is  given  to  within  about  25^  by  the  linear  analysis  at  the  fast 

shock  singularity  (3-28).   In  particular,  to  lowest  order  in 

2     2  2 

m  =  M.  -  1,  the  peak  magnetic  field  is  B^  =  sin  0  +  m  cosec  9 

and  the  period  27r/m.   But  note  that  these  oscillations  do  not 

strictly  correspond  to  any  steady  flow  oscillations  for  the 

entropy  is  still  zero  in  these  cases.   All  the  solutions  in 

these  two  limiting  cases  are  typical  of  the  "wave  train"  type. 

They  have  been  adequately  described  in  Section  3  and  examples 

are  given  in  Figures  6  and  7- 

There  is  just  one  feature  of  the  time  dependent  solutions. 
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Figure  10.  The  Cnief  Parameters  (0,€,M  )  on  Which  Solutions  Depend, 


showing  the  effect  of  the  artificial  viscosity,  that  should 
be  noted.   Before  breaking  occurs  there  is  no  need  of  artificial 
viscosity.   At  the  limiting  Mach  number  for  which  the  wave  just 
breaks,  one  obtains  solutions  like  that  on  the  left  of  Figure 
11  in  which  V  <=^  0  at  the  peaks  of  the  oscillations.   But 
introducing  viscosity  produces  the  solution  on  the  right. 
This  dramatic  change  is  understood  by  reference  to  the  phase 
diagrams  of  Figure  5(a).   The  region  of  complex  V-roots  leaves 
a  large  gap  in  the  possible  periodic  solutions,  which  can  only 
be  bridged  by  a  gain  of  entropy.   With  no  viscosity  the  solu- 
tion attempts  to  follow  the  solitary  wave  round  the  outside  of 
this  region,  but  when  viscosity  is  added  it  can  switch  to  the 
inner  solution  via  the  other  sheet,  as  it  must  after  breaking 
occurs . 

As  one  moves  away  from  these  special  cases,  predictions 
based  on  the  steady  state  theory  become  more  tentative  and 
one  has  to  rely  increasingly  on  the  numerical  solutions  of  the 
time  dependent  problem.   Let  us  consider  first  cases  in  which 
M„  <  Mp(0).   The  first  noticeable  change  is  the  introduction 
of  a  y-component  of  the  magnetic  field,  still  adequately  pre- 
dicted by  equations  (3-39)  and  (^.■^O).   Then  the  precursor 
begins  to  make  its  appearance.   Figure  12  shows  solutions  for 

fixed  values  of  9   and  M.  and  after  a  fixed  time,  but  with 

2 
four  different  values  of  e  .   The  behaviour  is  typical.   The 

oscillations  of  the  wave  train  have  periods  roughly  propor- 
tional to  the  mean  electron-ion  gyro-radius  but  amplitudes 
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2 
which  decrease  rapidly  with  decreasing  e  .   On  the  other  hand, 

the  precursor  scales  like  the  ion  gyro-radius  and  attains  Its 

2 
final  amplitude  within  a  factor  of  ten  change  In  e  .   From 

2 
there  to  e  =0,  shown  In  Figure  1J> ,    It  merely  Increases  in 

length. 

For  very  weak  waves  the  point  of  changeover  is  at 

tan  9  =   e~    -   e.      But  for  larger  M„  the  precursor's  appearance 

2 
is  delayed  to  smaller  e   (or  0) .      This  is  because  there  is  a 

wide  band  in  which  the  characteristic  roots  at  the  initial 

singularity  are  complex,  the  width  of  the  band  increasing  with 

M„  as  shown  in  Figure  3-   Thus,  for  comparison  with  Figure  12, 

when  9   =   tt/^  and  M.  =  1.^9  this  band  is  roughly  0.0^  <  e  <  O.56 

while  tan  6   =   e~^    -  e  at  e^  =  0.209-   Similarly  with  e   fixed 

at  0.275  X  10   ,  the  mass  ratio  for  a  deuterium  plasma,  the 

theoretical  angle  at  which  the  precursor  should  appear  when 

But  numerical  calculations  show  the 

following  precursor  lengths  at  t  =  10.0  (the  values  given  are 

the  ratio  of  the  precursor  length  to  the  length  of  the  whole 

solution  from  the  plasma -vacuum  Interface): 


M„  is  near  unity  is 
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It  must  not  be  forgotten  that  the  solutions  in  Figure  12 
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3. Or 

2.0 
B, 


60  80  100 

no  artificial  viscosity 


60  80  100 

artificial  viscosity 


Figure  II.    Solutions  Vi/hich    Just   Break,    Showing   the  Effect  of 
Adding   Artificial   Viscosity:  6^=  60*»,€2  =  |,  m^=  1.93.    The 


Time  is  50.0  . 
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2  4  6  8  10  12  14 

Figurei2.     Change  of  Field  Prof ile  With  Mass  Ratio:  0  =  60°,  M^=  I.  49.     The   ion   Gyro- 
Radius  is  the    Unit  of  Length,  the  Time  is  5.0. 
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are  very  time  dependent  and  that  the  time  at  which  they  are 
given  is  quite  early  in  the  development  of  the  precursor. 
Figure  1^  shows  how  such  solutions  change  with  time.   The 

precursor  becomes  very  irregular,  as  expected  from  the  almost 

2 
periodicity  of  solutions  near  the  initial  point.   As  e   is 

allowed  to  approach  zero,  the  higher  frequency  tends  to 

infinity  and  the  solutions  become  like  those  in  Figure  1^ , 

where  the  time  dependence  is  limited  to  a  gradual  extension  of 

the  precursor  ahead  of  the  main  wave  front.   These  were  obtained 

2 
by  setting  e   =  0  and  confirm  our  previous  assumption  that 

solutions  for  this  case  are  actually  the  limits  of  solutions 

o 
of  cases  in  which  e  — >  0  (except  for  9   =  Tr/2 )  . 

The  results  shown  in  Figure  13  merit  further  consideration. 
For,  since  none  of  the  solutions  of  the  steady  state  equations 
satisfy  the  initial  conditions  in  this  case,  it  might  be  thought 
that  they  have  little  quantitative  significance.   However,  the 
period  of  the  precursor  oscillations  is  given  by  equation  (J). ^3) 
as  30  compared  with  the  observed  first  period  of  2.9.   Also 
equation  {J>A^)    yields  for  the  ratio  (B^  -  sin  9):  Bp   the 
value  1.32  compared  with  the  observed  1.2  :  0.9-   Further  the 
dispersion  relation  of  (2.25)  yields  a  velocity  of  1.4  for 
oscillations  of  this  period  in  good  agreement  with  M„  =  1.42. 
For  the  oscillations  at  the  head  of  the  precursor  which  have 
a  period  of  about  1.2  the  dispersion  relation  gives  a  velocity 
of  2.8  in  good  agreement  with  the  observed  growth  of  the  pre- 
cursor. 
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Comparing  the  decay  of  the  precursor  (and  also  of  the  wave 
train)  away  from  the  main  wave  front  with  the  steady  flows  when 
friction  is  present,  tempts  one  to  say  that  the  time  dependence 
acts  like  a  friction.   That  there  Is  some  truth  in  this  can  be 
seen  as  follows.   Both  sets  of  oscillations  are  such  that  at 
any  point  a  fixed  distance  from  the  wave  front  the  oscillation 
grows  with  time.   Thus  transforming  any  variable  w  to  the  moving 
frame , 

w(x,t)  =  w(x  -  M„t,  t)  =  w(x,t), 

we  obtain  extra  terms  S.w  in  the  steady  state  equation,  for  which 
w~  B,w  =  a  >  0.   In  particular,  the  current  equation  (2.5)  gains 
a  term  aj,  which  constitutes  a  friction  when  a  is  constant,  as 
remarked  at  the  end  of  ^  ^•5- 

At  Mach  numbers  higher  than  Mp(0),  all  that  is  different 
is  the  form,  of  solutions  behind  the  wave  front.   At  one 
extreme,  e  =  1  or  0  =  7r/2,  the  solution  decays  monotonically 
to  the  final  constant  state  giving  the  discontinuous  steady 
state  solution  already  described  and  shown  in  Figure  7-   At  the 
other,  e  =  0  or  0  =  0,  the  solution  is  oscillatory  behind  as 
well  as  in  front  of  the  front.   The  period  of  these  oscilla- 
tions Is  small  near  Vl.pXP)    and  Increases  as  iyi„  increases,  as 
indicated  by  equation  (3.3^) •   Between  these  two  extremes  there 
is  a  wide  band  of  cases  in  which  there  are  complex  character- 
istic roots  at  the  final  singularity  and  the  solution  usually 
decays  rapidly  to  the  final  state.   Since  the  final  singularity 
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FigureK.  Development  of  Field  Profile  With  Time  for  on  Intermediate  Case: 

6=  75»,  €^=  0.01,  M^=  1.34.    The  ion  Gyro-Radius  is  the  Unit  of  Length  . 


In  all  these  cases  lies  on  the  shock  sheet,  there  is  always 
a  discontinuity  in  the  specific  volume  so  that  artificial 
viscosity  is  required  to  obtain  an  accurate  numerical  solution. 
The  accompanying  entropy  gain  means  that  there  is  a  non-zero 
fluid  pressure  behind  the  front  and  the  situation  at  the 
plasma -vacuum  interface  becomes  more  complicated.   To  main- 
tain overall  pressure  balance,  the  magnetic  field  imposed  in 
the  vacuum  to  obtain  these  solutions  has  to  be  considerably 
higher  than  that  at  the  fast  shock  singularity.   But,  as  we 
have  remarked  before,  a  contact  discontinuity  cannot  be 
maintained  when  9   <  7r/2  and  the  required  adjustment  takes 
place  through  a  rarefaction  wave.   The  situation  is  very 
similar  to  that  in  the  one-fluid  model  which  has  been  described 
in  detail  by  Bazer  [5]- 


Appendix  A.   Some  bounds  for  M»  and  proofs  of  two  statements 
in  the  text . 


For  simplicity  we  take  y  =  2,  p  =  0  and  denote  sin  6, 
cos  9,    and  B-^  by  s,  c,    and  B  respectively.   We  first  obtain 
bounds  for  M.  as  a  function  of  B  at  the  fast  shock  singular 
point.   Then  these  are  used  to  prove  statements  made  in  the 

text . 

2 
A  lower  bound  for  M„ : 

(A.l)      M^  >  c^  +  I  B(B+s) 

To  prove  this  we  need  to  consider  two  cases. 

(i)   0  <  rj  =  B/s  <  3:   Then  the  positive  sign  has  to  be  taken 


2 
in  the  expression  (3-19)  for  M„  so  that 


But 


M^  -  C^  =  T]  (l  -2c^+  [l-c^s^{Ti-lf]^   /3-11 


[l-c2s2(ri-l)2]2  >  i_l  (n-i)2  ..  1  (3-ri)(n+l 


so  that  also 

1  >  ^  (3-ti)(ti+1)  . 


Hence 


M^  -  c^  >  J  n(n+i)(i-2c2+  I)  -  I  b(b+s: 


(li)   11  >  3:   This  occurs  only  for  6   <  tt/^  and  the  positive 

sign  in  (5. 19)  gives  the  smaller  value  of  M. .   We  put  x  =  B-3s  >  0 
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and  expand  the  root : 

xM^  =  3s(2c^-l)(+x(3c^-l)-(x+3s)[l-c^(2s+x)^]2 

=  5s(2c^-l)+x(5c^-l) 

-  (x+5s) (l-4c^s^ )*[ i-ic^ (l-4c^s^ ) "^ (4sx+x^ )+•••] 

=  c^x+  (x+3s)(2c^-l)[ic^(2c^-l)~^(4sx+x^ )+...] 

2  2  —      2 
since  (1-4g  s  )^  =  2c  -1.   Clearly  all  the  terms  In  the  expan- 

2       2 
slon  are  positive  and  as  c   >  (2c  -l)  we  have 

'  .  .;      M^  >  c^  +  B(2s  +.ix)  =  c^  +  iB(B+s)  . 


2 
An  upper  bound  for  M.  when  the  singular  point  Is  on  the  Initial 

sheet : 


(A. 2)       s(M^  -  c^)  <  B^  . 

(1)      l<Tl<n      =1   +T^-      then   from    (5-19) 

s(M^    -    c^)    <   B(5-n)"^(l-2c^+l)    =   2Eyr)^(3-Tl) 
But    2   <    ri^(3-T])    If    (ti-i)(t]-ti^)(ti+ti^-2)    <    0    . 


(ii)   T]  >  "Hp,^  sin  6   >  1/3:   On  the  Initial  sheet  we  have 

2FV  +  G  >  0 
i.e.,      (6V-4)M^  +  2(B^-s^)  >  0 


-  91  - 


Substituting  for  V  from  (3- 17)  we  get 


;2B-3s)(M^-c^)  <  B(B^+l-2s^ 


As   2B-3S   >    (2/3^-1)3,    the    result    follows    if    l-2s^   <   2(/T-1)B^; 

i_ 
this  is  satisfied  for  all  n  >  n   if  (6+  kf^)^s   >  3s  >  1 . 


(ill)   n  >  ^  ^  sin  9   <  1/3:   It  is  easily  verified  that  if 


0  <  a  <  1, 


(y-a)-l[(l-a2)2_(i_y2)2j  <  ( i-a^ ) "^ (y+a) ,   0<y<l 


Thus  for  the  whole  range  ofB,  s<B<s+c~  ,  putting 
y  =  c(B-s),  a  =  2sc  into  (3.19)  and  using  the  above  gives 


(A. 3)     M^  -  c^  <  c^(l-2s^)  ^B(B+s) 


In  the  desired  range  we  therefore  get 


sCM^-c'^)  <  f(ri;^  +  %)B^   <  B' 


Statement  1  (in  g  3.2(b)).   At  the  fast  shock  singular  point 
on  the  initial  sheet  i-l-v  <  0. 

Since  2FV  +  G  >  0,  we  need  to  show 

(^x-v)  (2FV+G)  =  [6m^V-2(2m|+s^-B^)]  (V-c^M^^)  -  2B^V 


-  -2B  ^{v[B^-s(M^-c^)]  + 
+2s(l-c^M^^)(l-sB"^)[M^-c^-iB(B+s)](<  0, 


where  we  have  used  (3.I7)  to  substitute  for  V. 
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But  this  follows  Immediately  from  (A.l)  and  (A. 2)  above 


Statement  2  (in  ^3-3(c)).   The  quadratic  (3.13)  yields  real 
roots  for  V  on  the  B^-axis  between  the  origin  and  the  fast 

shock  singular  point . 

2 
The  discriminant  of  (3-13)^  G  -FH,  is  a  quartic  P(B^) 

in  B^  of  the  following  form 

(A. 4)  p(x)   =    (x^    -   B^f  +   3[(cx  +x^f-   b^] 

p  op  pii?P?  ? 

where  a     =  21^^  +    s      ,      b      =  M^  +   s    (2M^   -   c    )   +  x      , 

and  X     =   s(M^    -   c^)   >   0. 

o  ^    A  ^   — 

We  note  first  that  P(0)  >  0;  this  requires 

p    pp      h  '?   o  '^   "? 

(2M^  +  ^    )      >   3(M^  +  2s^M^  -  s^c^) 

ii     7^  o         h  "^  '^  p    pp     pp 

i.e.,      M^  -  2s  M^  +  s   +  3s  c   =  (M^  -  s  )   +  3s  c   >  0  . 

We  also  know  that  P(sin  0)  >  0  and  P(B^)  >  0.   Further,  if 
P(^)  <  0  for  ^  >  0  then  P(-^)  <  P(^)  <  0.   Thus  there  can  be 
at  most  one  interval  on  the  positive  real  axis  in  which 
P(x}  <  0.   We  have  to  show  that  this  lies  beyond  B^.  It  exists 
only  when  P(x)  =  0  has  four  real  roots j  let  them  be  :: 
X-,  <  Xp  <  x^  <  ^ii  •   Then  we  have  only  to  show  that  x^,  >  B^  =  B 
if  x^  >  x^  >  0. 
From 

y       X   -  0  ,    ^_  XX   =  A„,   X  X  X  X.  =  A.  >  0 

i=l   ^         i<j   ^  J    "^    1  ^  p  ^    H 
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we  obtain 


But 


xjx^  +  x|x^(x^+  x^)  +  A^x^x^  -  A^  -  0 


3xj^  ^  x^  +  x^(x^+x^)  =  A^(x^Xi^)  ■'■-  Ag  ^  A^xj;^  -  A^ 


so  that 

4       2 
3Xi^  +  A^Xj^    -  A^  >  0  . 

Substituting  for  Ap  and  A|.  from  (A. 4)  we  find  the  roots  of 
this  quadratic: 

^  (4m^+  5s^  -  3  t  [28m^  -8(1+2c^)M^  +  s^  +6s2+  9]*  I 

2 
Since  M»  >  1,  the  expression  under  the  radical  sign  Is  greater 


A 
than 


'A   2 

^4  —  --  ^ 


5M?  -  i  -  2c2)2 


2 
and  Xi|  cannot  be  less  than  the  positive  root.   Thus 


12x^  >  18m^  +  3  -  l4c^  >  9B^  +  9sB  +  3+  4c^ 


by  (A.l).   Also 

9sB  +  3  >  3B^  If  0  <  2B  <  3s  +  (9s^  +4)2  . 

This  Is  clearly  satisfied  for  9   >  ir/^,   when  B  <  3s;  and  when 
9   <  7r/4  we  have 

B  <  s  +c"^  =  s  +  (1  +s^  +  s^+.  .  .  )^  <  s  +  (l+2s^)2 

so  again  It  Is  satisfied.   Thus  we  have  the  desired  result 

X4  >  B. 
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Appendix  B-   Energy  Identities  and  stability  conditions  for 
the  difference  equations. 

We  shall  consider  here  only  those  cases  in  which  no 
breaking  occurs.   Thus  a  smooth  solution  to  the  differential 
equations  can  be  assumed  and  no  artificial  viscosity  is  needed. 

It  Is  assumed  that  a  solution  is  required  for  a  finite 
time  Interval  0  <  t  <  T,  and  that  there  is  a  finite  number  R 
such  that  all  solutions  are  constant  for  x>R,  0<t<T. 
This  removes  the  spurious  difficulty  of  unbounded  energy 
Integrals.   Then  we  define 

N 
(B.l)      (u,v)  .  =  h  ^_  u  V   ,   ||u  IL  =  (u,u)  . 

for  real  mesh  functions  u,v,  where  N  =  [R/h]+2.   When  H   =   0 
the  subscript  will  usually  be  dropped.   The  manipulations  mainly 
Involve  summation  by  parts  for  which  one  needs  the  formulae 


:b.2 


:u,D_v)_g  =  -(V^)i-l  +  "n+1^N  -  "^-l^^-l 


V^i  -  -\^.^'^u+-i  ^  ^N^N+1  -  "ri  ' 


where  D  u.  =  h~  (u  .  -,  -  u.)  and  D_u  .  =  h   (u  .  -"._-,).   Finally, 
we  denote  by  c  ,  c.,  and  c  the  adlabatlc  sound  speed,  Alfven 
speed,  and  magneto-sonic  speed  on  the  Lagranglan  coordinate 
space,  i.e., 

(B.3)      c^  .  YA/V^+^  ,  ol   =   B.B/V  ,   c^  =  c^  +  c2  . 

-  95  - 


1 .   A  non-llnear  Identity  for  the  transverse  case. 

The  difference  equations  In  this  case  are,  for  j  =  0,1,..., 


1.5)     V  .    =  V  .  +  kD  ,  u  . 


(B.6)      V"+^  B^+1  -  D  D  B^+1  =.  1 


(B.7)      a'}^^  =  a" 


We  define  the  energy  as 


;b.8)  e  =  I  ||u  1|2  +  |||d^b||2  +  |||v^B||^  +   (t-i)-^p,v) 


Then  taking  the  difference  of  (B.6)  at  n+1  and  n,  multiplying  by 

n+1    n 
B.   +  B.,  and  summing  over  j,  we  obtain 

J  J 

(v"+V+^b"+Vb")-(v"b",b^+Vb")  =  (d^d_[b"+1-b"],b^+V) 
=   -(dJb^+1-b"],d_[b^+Vb"])^-(b^+i+b^)d_[b^+1-b^] 

+  + 


If  the  homogeneous  boundary  condition  B=0orSB=0  1s  Imposed  at 
X  =  0.   Multiplying  (B.5)  by  B  .B .  .  and  summing  gives 

(^n+l_  ^n^  gn^n+l^  ^  k(D^u"+\  b"b"+1) 
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which,  when  subtracted  from  the  equation  Immediately  above, 

yields 

(B.9)    I  KV^-^^J^B^+I  ||2^||j^^3n+l  ||2_  j,  ^^n^i^n  |,2_  y^^^n  |,2 

=  -k(D^u^+\  b"b^+1)  . 

Now,    since   A.    Is   Independent    of  n, 

/      tn-1/    n+l-.n+l        '^n,,nN  r^-n/^n+l     ^nx 

(y-1  )       (P  .     V  .         -  p  .V  .  )    =    -p  .  ( V  .      -  V  .  ) 

=    -kp  .    D^u. 

where  p^  Is  a  pressure  Intermediate  between  p .  and  p .   .   So 

J  J  J 


/^  -,  ^\     I      n  \  -If  /  n+1  -,^n+l  X  ,  n  -.nx  -,     ,  /_   n+1  ~n. 
(B.IO)     (y-1)   [ (p    ,V    )-(p  ,V  )]  =  -k(D  u    ,p  . 


Lastly,  from  (B.4)  after  multiplication  by  u .  -,  +  u  ._|_-|^ 

(B.ll)    ||u"+l  11^  -  ||u^^  11^  =.  -k(u"+l  +  u^D_[i(B")^  +  p^])l 

^  k(D^[u"+l+u"],i(B^)2  +  p") 


If  the  boundary  condition  u  =  0  at  x  =  0  Is  assumed  and  u  =  0 
Initially. 

Thus  If  the  energy  Is  modified  by  putting 


,.12)     L  =  E  +  ik(D^u,  \-^   +  p) 


we  obtain  the  following  relation  from  (B.9)-(E.11) 
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(B.13)     L     -  L  =  2k(D_^u    ,  2LB    -  B  J   +  [p   +p  -2p  ]) 

We  shall  not  carry  the  consideration  of  this  identity  any 

further  except  to  note  that  the  right  hand  side  will  usually 

be  0(k  )  and  to  see  under  what  conditions  L  is  positive  definite 


I.e.  , 


ik|(D_^u,iB^+p)|  <  E 


Now     _  ^k  I  (D^u,  iB^+p)|  <  (k/h)  1 1  ^B^  +  p  1 1  ||u 


2  2    2 

and  the  inequality  (ax-by)   <  (a+b)(ax  +by  )  for  a,b  >  0  gives 


iB^+  p)^  <  I   V-^(b2+  2(y-1)p)(b2v  +  2(y-1)-V: 
=  ^  (c2+2(l-y-l)c2)(B2v+2(Y-l)-V; 


Thus  L  is  positive  definite  if 


;b.14)     (k/h)max(c^  +  2(1-7-^)0^)2  <  2 


2 .   A  Sufficient  condition  for  stability  and  convergence 

To  prove  the  convergence  of  the  difference  approximations 
to  the  assumed  solution  of  the  differential  equations  as 
h,k  — >  0,  it  is  necessary  to  linearise  about  this  solution. 
We  begin  with  the  transverse  case,  and  for  definiteness  assume 


that  the  set  of  functions  u ,  B,  V,  A  constitutes  a  smooth 
solution  of  the  differential  equations  with  B  and  u  prescribed 
at  X  =  0.   We  assume  also  that  V  >  V  >  0.   Further,  let  u,  B, 
V,  A  be  the  solution  of  the  difference  equations  (B.4)-(b.7) 
for  the  same  boundary  and  Initial  data,  and  put  u  =  u-u,  etc. 
Then  u=B=Oatx=0  and  A  =  0  everywhere.   Also 

~2     2     ^        ^2 
B   -  B  -  2BB  +  B   , 


and 


V  ^  -  V  '^  =  -yV  ^^■^^•^  V  +  O(V^)   If   |v|  <  iV^,  say 


Hence,  If  we  denote  by  T-,  ,  Tp ,  T^  the  truncation  errors  obtained 
by  substituting  the  solution  u,  B,  V,  A  Into  the  difference 
equation  (B.4)-(b.6),  we  get 

(B.15)     Uj.+i  =  u^.^^-kD_[B.^^B.^3_-(c3^j.^^)  V.^^]  +  T^ 

+o[(v"+b")2] 

(B.16)     V^+1  =  V^  +  kD  G^+1  +  Tp 

^^     JJ      JJ      +-J      JJ      3 
This  time  we  take  for  the  energy  E  and  Its  modification  L, 


E  -  i  iiu  ir  +  i  iiD  B  ir  +  i  iiv^B  II   +  i  iic  vir 


(B.18) 


L  =   E  +   ik(D,u,    BB    -    c   v; 
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It  Is  also  convenient  to  Introduce  the  quantity 

g  =  lul  +  V^ IbI  +  c  Ivl 
for  which 

(B.19)    hg^  <  ||g  11^  <  6e  . 

Then  in  a  manner  similar  to  that  above  we  obtain 

(B.20)  ||u  ||-||u||       =k(DLU         +u],BB-(c)V)+R. 

with  R^   =   OHT^   +    (g    )    ,    g    +  g 

(B.21)  ||cflv"+l||2-||c^v^||2  =  k(D^G"+\(c^+M2-n+l^^^n^2^n)^^^ 

with  R^    ^    (T2,(cf  l)"v"-*-^+    (c2)'v")-([cf  l)2-(c^)2]v^v"+l) 

=  0|(T2,  g"  +  g"+^)  +  h(i^  i"+i) 

(B.22)  ||(V"+^)*   b"+^||2+||d   S^+l||2-    ||(V")*B^||2-||D  g"||2 

+  + 

=  -k(D^u"+\  b"+^b"+^  +  bV)  +  R^ 

with     R3  =  -(v"+i-  v",  b"b"+1)-  (b"+1-  b^,  v"+V+v"S^+^) 

=  0^(T2+T-^+(g)+(g    ),g+g    )+h(g,g    )| 
Collecting  these  results  yields 
(B.23)    l"^^  -  l"  =  R^  +  R^  +  R-^  . 
It  is  now  necessary  to  ensure  that  L  is  positive  definite  and 
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then  show  that  It  tends  to  zero  with  h.   For  the  first,  we  have 

ik| (D^u,BB-c^V) I  <  (k/h) I (u,BB-c^V) I 
n  II  ^  II  2    1  /  /  n2  1 1  ■^   2^^  II  2 

<  i  iiu  ir  +  i(Vh)  iiBB-c^v  ir  • 

2  2      2 

But  from  the  Inequality  (ax  -  by)   <  (a+b)(ax  +  by  ) ,  we  get 

(BB  -  c J)^  <  (bVv  +  cl){VB^   +    c^V^) 

=  c^(VB^  +  c^V^)  . 

Thus  L  Is  positive  definite  if 


(B.24)     (k/h)  max   c  <  1 
(x,t) 


In  order  to  estimate  the  right-hand  side  of  (B.23),  we 
observe  first  that,  if  the  truncation  errors  T-,  ,  Tp ,  T^  are 
bounded  by  T^,  then 

{B.25)    S"+^  =  0(T„  +  r   +  (S")2)  . 

So  far  as  u  and  V  are  concerned,  this  follows  directly  from 
the  difference  equations  (B.15)  and  (B.lb).   And  for  B  we 
rewrite  (B.17)  in  the  form 

'^n+l    /~n+l      v-^n+l         n+1  ^n+1 
MB     =  (V.   -  D  D  )B.    =  T-^  -  B.    V.    , 

B   =  B,.  =  0 
o     N 

and  notice  that  the  operator  M  satisfies  ||mI|  >  ^V  when 
|V    I  <  iv  .   Thus  we  conclude  that  there  exists  a  positive 
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constant   a.   Independent  of  h,  such  that 


(B.26)    T^+  g"<a   Implies  Iv"""^^!  <  ^V^ 

'^n+1 
and  hence  J  from  the  estimate  for  g    ,  that 


-^"  =  o[tm  +  (^h+  C^+f 


For  our  difference  scheme,  T  =  0(h),  and  If  the  stability 
condition  Is  satisfied,  E  =  0(L  ).   So  there  exists  a  constant 
K,  Independent  of  h,  such  that  If  (B.26)  Is  satisfied 

(B.27)     k-l(£"+l-  £^)  <  K|h5+h2(£^)*+L'%h-5/2(£^)^/2J 

Now  let  f(t)  be  the  solution  of 

(B.28)    ^  =  2K(h  +  fi  +  f),    f(0)  -  0  . 

Then  f  Is  bounded  by  some  constant  F  for  0  <  t  <  T,  and  we 
can  show  by  Induction  that 

(B.29)    l"  <  h  f (nk)  ,     0  <  nk  <  T  . 

For,  from  the  assumption  of  (B.26)  and  (B.29)  for  n  =  0,1,..., m. 
It  follows  that,  with  f  denoting  f(mk), 

,  -l/:^m+l   ^ms  ^  ^^,  4/,   ^i   ^   ,  i  ^3/2s 
k   (L    -  L  )  <  Kh  (h+ f2  +  f  +h2  f-^/ "") 

—  ^     m   m      m   ^ 

<  h  (df/dt)      <  h^k"-^(f  ^,  -  f  ) 

-  '       ^t=mk  -        m+1    m^ 

-1         '^m+l     k 
for  all  h  <  F   .   Thus  L    <  h  f^^^  which  Implies,  by  (B.19), 
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that  g"^"^-^  =  0(h-^/^).   So  for  s^jf flclently  small  h,  (B.26) 
can  be  satisfied  for  n  =  m+1  and  the  Induction  continued. 
The  final  result  follows  If  L  =  0  Initially  (in  fact,  one  may 
have  L  =  0(h)). 

The  extension  of  this  result  to  the  case  of  arbitrary 
field  orientation  Is  quite  straightforward.   The  mass  ratio 

has  no  effect  and  the  only  change  Is  that  the  Alfven  speed  is 

2      2     2       2 
now  given  by  c.  =  (Bp  +  B^  +  cos  9)/V. 
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